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Abstract 


This research focused on incorporating stability and control into a multidisciplinary de- 
sign optimization on a Boeing 737-class advanced concept called the D8.2b. A new method 
of evaluating the aircraft handling performance using quantitative evaluation of the sys- 
tem to disturbances, including perturbations, continuous turbulence, and discrete gusts, is 
presented. 

A multidisciplinary design optimization was performed using the D8.2b transport air- 
craft concept. The configuration was optimized for minimum fuel burn using a design range 
of 3,000 nautical miles. Optimization cases were run using fixed tail volume coefficients, 
static trim constraints, and static trim and dynamic response constraints. A Cessna 182T 
model was used to test the various dynamic analysis components, ensuring the analysis was 
behaving as expected. Results of the optimizations show that including stability and con- 
trol in the design process drastically alters the optimal design, indicating that stability and 
control should be included in conceptual design to avoid system level penalties later in the 
design process. 
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Chapter 1 
Introduction 


Current trends in transport aircraft conceptual design search for technologies that pro- 
vide incremental performance improvements, many focusing on engine efficiency and drag 
reduction technologies resulting in reduced mission fuel burn. The National Aeronautics 
and Space Administration (NASA) is pushing for future transport concepts that go beyond 
incremental reductions in fuel burn, noise, and emissions. Achieving the aggressive reduction 
goals set forth by NASA will require a complex, multidisciplinary design process, integrating 
each discipline early in the design phase. Historically, aircraft conceptual design has included 
aerodynamics, sizing, weight estimation, propulsion, and mission analysis. Stability and con- 
trol, and structural design, among others, often have been left with a fixed design that is 
frequently a sub-optimum solution for those disciplines, leading to system level penalties. To 
achieve the aggressive performance goals set by NASA’s Fixed Wing (soon to be Advanced 
Air Transport Technologies) Project, a full integration of all disciplines is required. 

Aircraft drag consists of a combination of lift- independent and lift-dependent, mostly 
induced when sub-transonic, drag. In commercial aircraft, the trend is to increase wing 
span to achieve a lift-dependent drag benefit. Increasing aircraft span, and aspect ratio for 
a fixed wing area, reduces lift-induced drag, but aspect ratio is currently limited both by 
structural (flutter) and dimensional constraints (airport operations). For subsonic flight, the 
total aircraft lift-independent (parasite) drag consists mostly of skin-friction and separated 
flow pressure drag [1]. With fairings and surface blending, much of the pressure drag can be 
minimized, leaving skin-friction drag as the largest contributor to parasite drag. Effectively, 
if the wetted area of the aircraft can be minimized, the skin-friction drag will be minimized, 
resulting in reduced parasite drag. Numerous designs have reduced this wetted area by using 
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tailless configurations, such as the Convair F-106, the Convair B-58, the Messerschmitt 163B, 
and the Northrop Grumman B-2 [2,3]. All of these configurations have low aspect ratios, 
resulting in reduced cruise efficiency and low-speed aerodynamic performance. It would be 
beneficial if both the induced drag and parasite drag could be reduced simultaneously. 

For a conventional configuration, eliminating lift produced by the empennage can reduce 
aircraft induced drag, even if the tail is producing positive lift. Aerodynamic efficiency is a 
function of span loading, and the horizontal stabilizer aspect ratio, for stall considerations 
related to safety, is less than that of the main wing. As induced drag coefficient is a function 
of aspect ratio, any lift produced by the empennage with span loading less than the main 
wing is less efficient, resulting in an induced drag penalty. 

Adjusting the aircraft static margin is a way to reduce induced drag on the empennage. 
As the static margin decreases, the aircraft will first become decreasingly statically stable, 
then neutrally stable, and eventually become unstable as the static margin becomes negative. 
Decreasing the static stability reduces the load on the horizontal stabilizer, allowing for 
greater aerodynamic efficiency and reduced structural weight. Many modern aircraft, both 
commercial and military, are designed with relaxed static stability, the intentional decrease 
of static margin, which requires some form of stability augmentation in order to reduce 
the trim-drag penalty [4-6]. As stated by Raymer, “[A] modern and sophisticated aft-tail 
aircraft is designed to a slight level of instability so that it normally flies with an upload, 
not a download on its tail. This is the very reason that computerized flight control systems 
with artificial stability were developed and put into production” [1], 

Classical aircraft design uses volume coefficients when sizing an aircraft’s empennage, 
which tends to produce conservative estimates for the stabilizer areas [7]. As a result, the 
surface area for the horizontal and vertical stabilizers exceed the necessary area for both static 
and dynamic stability. Additionally, this fails to take into account the benefits of augmented 
stability provided by the active control systems in modern aircraft. As mentioned previously, 
a majority of parasite drag results from skin friction over the wetted area of the aircraft. 
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Augmenting the aircraft’s stability will allow the stabilizing surfaces to be reduced in size, 
thus reducing the total wetted area of the aircraft. In doing so, not only can the induced 
drag be reduced from the relaxed static stability, but the total parasite drag can be reduced 
from the decrease in wetted area. 

The benefits of designing for relaxed static stability are not unbounded, however. With 
additional instability, the workload of the active control system increases to a point where 
it can no longer provide adequate dynamic stability to the aircraft. Although the total drag 
may continue to decrease, the design becomes infeasible and uncontrollable. Additionally, the 
active control system needs to stabilize the aircraft dynamics in such a way so as to provide 
acceptable handling qualities. For manned aircraft, MIL-F-8785C (Ref. [8]) was formulated 
from years of handling-qualities research as a guideline which has been used extensively 
in both military and commercial aircraft. These handling-quality guidelines were developed 
using flight tests and simulators which were influenced by pilot opinion. The current handling 
qualities specifications, MIL-STD-1797A [9] which evolved from MIL-F-8785C [8], require 
natural frequencies and damping specific to each system linear mode, a limitation if the 
typical linear modes do not exist. For unconventional, advanced configurations with stability 
augmentation systems, this can be the case. 

Multidisciplinary analysis is required to capture all the coupled effects between aircraft 
sizing, aerodynamics, weight estimation, propulsion, and stability and control. Presented 
here is a methodology for integrating stability and control into the conceptual design pro- 
cess, where emphasis was placed on using conceptual level design parameters and reduced 
inputs for the controller. Atmospheric disturbances and perturbations were used to stress 
the flight dynamic system with static and dynamic constraints placed on the system re- 
sponse, used to size the configuration in a design optimization. Using this methodology, 
the stabilizers were sized using static and dynamic constraints with augmented stability, a 
physics-based approach, instead of the empirical tail volume coefficient. Multidisciplinary 
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design optimization was used to minimize fuel burn on a D8 advanced transport aircraft 
concept, shown in Fig. l.l 1 , including stability and control in the design analysis. 



Figure 1.1: Conceptual sketch of the D8 advanced transport aircraft concept. 

The D8 geometry is a Boeing 737 class unconventional geometry currently being studied 
by the Massachusetts Institute of Technology (MIT) and NASA, with key features including a 
double-bubble fuselage and rear embedded engines using boundary layer ingestion (BLI) [10- 
13]. Two concepts were developed in Ref. [10], both designed for a 3,000 nautical mile mission 
focusing on fuel burn minimization. The first concept was a current technology concept 
showing the configuration benefits alone, and the second was an advanced N+3 (now plus 
three generations beyond current commercial aircraft 2 ) technology concept with expected 
entry into service near 2035. The current technology D8 concept was used in this research. 

The process chart of Fig. 1.2 shows a high level overview of the analysis methodology 
used to incorporate stability and control into the conceptual design process. A combina- 
tion of the NASA Flight Optimization System (FLOPS) [14] and custom scripts defined 
the geometry and are discussed in Sections 3.1, 5.2.2, and 5.3.1. The stability and control 
derivative and aerodynamic analyses used the vortex lattice code Athena Vortex Lattice 
(AVL) discussed in Section 5.1, and the flight dynamics and system disturbances were pro- 
gramed using MATLAB®, discussed in Chapter 4. The mission analysis was performed 

1 http : / / www . nasa . gov/ topics/aeronaut ics/ f eatures/f uture_airplanes_index . html 

2 http : //www . aeronautics .nasa. gov/nr a_awardees_ 10_06_08 .htm 
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using FLOPS, where the mission fuel burn, including reserve fuel equal to 5% of the total 
mission fuel, was calculated. Each analysis component of FLOPS is described in detail in 
Section 5.2. The analyses were integrated into a multidisciplinary analysis software called 
ModelCenter®, created by Phoenix Integration, 3 and the system was optimized with the ob- 
jective of reducing total system fuel burn. The multidisciplinary framework and optimization 
methodology are discussed in Section 5.3. 



Figure 1.2: High level process chart of the multidisciplinary analysis used in this research. 

It was desired to test different pieces of the multidisciplinary design analysis individ- 
ually. Stability derivative flight test data was available from Ref. [15] for a Cessna 182T 
aircraft, shown in Fig. 1.3. The stability derivative data was used as a validation case for 
the aerodynamic analysis, and then the model was used throughout the methodology devel- 
opment to exercise each analysis piece as a sanity check, verifying the flight dynamics and 
disturbances were producing the desired behavior. 

Five optimization cases were run with the current technology D8 geometry. The first 
case was optimized using traditional fixed tail volume coefficients. Static trim constraints 
were then added to the optimization and the fixed tail volume coefficients removed from the 
geometry definition. Case two had a fixed static margin, and case three was allowed to vary 

3 http : / / phoenix- int . com/ 
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Figure 1.3: Picture of a Cessna 182T from the Pilot’s Operating Handbook [16]. 

the static margin. Optimization cases four and five added dynamic performance constraints 
to the optimization, with case four having fixed static margin and case five having varying 
static margin. Results indicate that the dynamic response constraints drastically altered 
the optimal designs, actively constraining the design space, highlighting the importance and 
benefit of incorporating stability and control into the conceptual design process. 
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Chapter 2 

Review of Literature 

2.1 Flight Dynamics and Control in Conceptual Design 

Classical conceptual design typically focuses on the interaction between the disciplines 
of aerodynamics, sizing, weights estimation, propulsion, and performance [17]. Any inclusion 
of stabilizer or control surface design during conceptual design is often limited to estimating 
sizes from historical data, assuming control effectiveness is proportional to the area and 
moment arm [7,18]. Often times, flight dynamics and control (FDC) and handling qualities 
are examined after the aircraft geometry and structural properties have been defined, which 
leads, inevitably, to sub-optimal designs or even configurations with deficient flying qualities. 
FDC is incredibly important when it comes to the overall safety and certification of an 
aircraft, [7, 19] and deficient handling qualities lead to reduced aircraft performance, large 
cost increases, and delays as the configuration must be re-evaluated or redesigned. 

Collaboration between traditional conceptual design disciplines and flight dynamics and 
controls is essential when trying to properly size stabilizer surfaces for advanced concepts 
that lie ontside the stabilizer-sizing empirical databases, especially when exploring the design 
space with reduced static margin, or relaxed static stability (RSS). RSS is the intentional re- 
duction of static margin with the objective of obtaining performance gains. When designing 
with RSS, it is possible for aircraft performance to be increased through reduction of wetted 
area drag, trim drag, and tail weight [7,18,20]. For a transport aircraft with conventional 
stability margins, the horizontal tail accounts for 20-30% of the aircraft-lifting surface and 
approximately 2% of the aircraft empty weight [18]. Any reduction in size of the horizontal 
stabilizer can provide a significant benefit in reduced drag and aircraft empty weight. How- 
ever, any relaxation of the stability margin has a detrimental effect on the aircraft’s handling 
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qualities that must be addressed [18]. It is this correlation between performance gains and 
degraded handling qualities that make it essential, and potentially extremely beneficial, to 
incorporate flight dynamics and control into the conceptual design phase, especially when 
evaluating large and complex design spaces using an optimizer. 

Perez, Liu, and Behdinan in Refs. [7] and [18] incorporate FDC into a multidisciplinary 
design optimization (MDO) scheme where a manned transport aircraft was allowed to be 
designed with RSS, taking into account handling qualities as specified in MIL-F-8785C [8]. 
From this optimization scheme, significant geometry changes occurred compared to using 
the traditional design approach, where aircraft performance was only a function of the aero- 
dynamic forces. Specifically, the optimizer recognized a benefit from RSS and reduced the 
static margin, along with moving the wing apex location and, thus, the center of gravity 
(CG) location by reducing the horizontal stabilizer area (down 28%) and control surface 
area [7,18]. By integrating stability and control into the design optimization loop, adequate 
handling qualities were ensured with feedback control augmenting the stability. Additionally, 
reduced control deflections were necessary for trim [17]. 

References [17] and [18] examined the longitudinal dynamics only, a limitation in the 
ability to fully design the system stabilizers and control surfaces. The lack of fully-coupled 
aircraft dynamics does not allow the vertical tail design to be incorporated in the MDO. As 
stated previously, sizing of the vertical tail is traditionally calculated using a vertical tail 
volume coefficient, developed from historical data. The dutch roll mode, which is primarily 
damped with the vertical tail, is made sufficiently stable simply by the vertical tail area and 
moment arm when using the volume coefficient [4] . With stability augmentation, this sizing 
will be less than optimal and potential performance gains are lost when the lateral/directional 
modes are not considered as part of the MDO. In 2006, Perez et al. expanded upon their 
research and incorporated both lateral/directional and longitudinal dynamics into the MDO, 
but the stability augmentation system only used Single-Input, Single-Output (SISO) feed- 
back controls with simple approximations of the common aircraft modes [7]. This does not 



account for an aircraft’s potentially complex, coupled dynamics, as would be the case with 
some of the next generation, advanced concepts being explored by the National Aeronautics 
and Space Administration (NASA). 1 Additionally, the gain direction and value was chosen 
to guarantee closed-loop stability but gave no indication of an optimal solution. Although 
able to produce results that showed a performance gain similar to their previous work, in- 
corporating MIMO feedback control gives the potential for more optimal feedback gains on 
a general configuration for all flight modes, and thus an optimal solution providing greater 
performance gains. 

In a completely different technique, Morris et al. used the method of linear matrix 
inequalities (LMI) to place constraints on the maximum actuator deflection, actuator rate, 
and pole placement limitations [21]. This method relies heavily upon the work of Boyd [22] 
and Kaminer [23] to place constraints on the static feedback gain matrix, K, to obtain desired 
handling qualities. Morris expanded his work in Ref. [24] by translating the MIL-STD-1797A 
guidelines into state variance constraints to be used in the development of a state feedback 
control law using optimal control. A very unique solution to the problem of incorporating 
FDC into the conceptual design process, the methodology was still linked to a standard set 
of linear dynamic modes and their associated natural frequency and damping. 

References [25-28] describe the SirnSAC project using the CEASIOM software which 
took a different approach than the methods described by Perez et al. and Morris et al. The 
SimSAC project used higher order tools, such as computational fluid dynamics (CFD), to 
iterate a conceptual design. The project had good success showing some benefits of relaxed 
static stability, but the higher order tools reduce the ability to explore a large design space 
with numerous varying geometric parameters. The optimization time using CEASIOM was 
on the order of weeks instead of the much faster methods described by Perez and Morris, 
and the methods presented in this research. 

^ttp : / / www . nasa . gov/topi cs/aeronaut ics/f eatures/f uture_airplanes_index . html 
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2.2 Handling Qualities and Stability Guidelines 


For years, research has been conducted to obtain a better understanding of the handling 
qualities of any particular aircraft, taking into account the human element, so a standard set 
of requirements could be formed. The earliest handling qualities requirement was published 
in 1965 by Westbrook and simply stated, “During this trial flight of one hour it (the air- 
plane) must be steered in all directions without difficulty and at all times be under perfect 
control and equilibrium” [29]. Although simple in nature, this less than clear statement has 
over time evolved into a much more extensive quantitative and sophisticated set of criteria 
that constantly change with every new class of vehicle [29]. The evolution was far from 
simple, however. Extensive tests in both simulators and in flight were used to quantify the 
human response which would allow for control system design to provide specified handling 
qualities [30,31]. 

It has been quite difficult to establish a standard set of guidelines for control systems 
that encompasses all types of aircraft and the various nuances of each individual pilot. The 
pilot’s opinion of how the aircraft handles is inevitably skewed by more than just how the 
aircraft handles, but more how it “feels.” “He or she will be influenced by the ergonomic 
design of the cockpit controls, the visibility from the cockpit, the weather conditions, the 
mission requirements, and physical and emotional factors” [4], A systematic method of 
quantifying test results was required to catalog all the data. The most widely accepted 
standard description for pilot rating handling qualities that establishes a quantitative scale 
was developed by G. E. Cooper and R. P. Harper, Jr. in 1969, and is called the Cooper- 
Harper scale [4,5,29-32], Shown in Table 2.1, the Cooper- Harper scale provides a way to 
correlate pilot opinion ratings with the aircraft dynamic model and determine some analytical 
specifications for good handling qualities [4], 

Cooper and Harper grouped the pilot rating scale into three separate levels that would 
describe the dynamic and control characteristics of an aircraft, described in Table 2.2. These 
guidelines can be applied to a specific flight mode where the three flying levels are correlated 
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Table 2.1: Cooper-Harper scale [4], 


Aircraft 

Characteristics 


Demands on Pilot in Selected Pilot 

Task or Required Operation Rating 


Flying 

Qualities 

Level 


Excellent; highly 

Pilot compensation not a factor 

1 


desirable 

for desired performance 

1 


Good; negligible 
deficiencies 

as above 

2 

1 

Fair; some mildly 
unpleasant 
deficiencies 

Minimal pilot compensation 
required for desired performance 

3 


Minor but annoying 
deficiencies 

Desired performance requires 
moderate pilot compensation 

4 


Moderately 

objectionable 

deficiencies 

Adequate performance requires 
considerable pilot compensation 

5 

2 

Very objectionable 
Vint, tnlprpmlp 

Adequate performance requires 

6 


ULIU uvJlUl CIjUIU 

deficiencies 

extensive pilot compensation 


Major deficiencies 

Adequate performance not 
attainable with maximum 

7 


tolerable pilot compensation; 
controllability not in question 




Major deficiencies 

Considerable pilot compensation 
required for control 

8 

3 

Major deficiencies 

Intense pilot compensation 
required to retain control 

9 


Major deficiencies 

Control will be lost during some 
portion of required operation 

10 
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Table 2.2: Flying quality levels [31]. 


Level Description 

1 Flying qualities clearly adequate for the mission flight phase 

2 Flying qualities adequate to accomplish the mission flight phase but with some 
increase in pilot workload and/or degradation in mission effectiveness or both 

3 Flying qualities such that the airplane can be controlled safely but pilot workload 
is excessive and/or mission effectiveness is inadequate or both 


to the natural frequency and damping ratio as shown in Fig. 2.1. Using these ratings, spec- 
ification MIL-F-8785 was written to provide dynamic stability requirements, both natural 
frequencies and damping coefficients, used both by military and civilian aircraft today. All 
of the qualities described have been related to how the aircraft handles according to the 
pilot’s opinion of the workload and the degradation of mission effectiveness. 



0.1 0.2 0.4 0.6 0.8 1.0 2.0 

Short-period damping ratio, £ 

Figure 2.1: Short-period flying qualities [31]. 
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Chapter 3 

Configuration Description and Mission Summary 


Two geometries were used in this research and are described in this chapter. The 
D8.2b, a current technology, span-restricted, twin-engine variant of the advanced transport 
aircraft described in Ref. [10], was used in a multidisciplinary design optimization described 
in Section 5.3. Throughout the methodology development and integration, a Cessna 182T 
model was used to test each component ensuring each was performing as expected. The 
Cessna 182T geometry is described in Section 3.2. 

Considerable effort was used to ensure the accuracy of the geometries used in this re- 
search with the original geometries. The geometries are described in the following sections 
along with the original geometry sources of the information, including drawings when appli- 
cable. 

3.1 D8.2b Geometry and Mission Definition 

The D8 series advanced concept tube-and-wing transport aircraft described in Ref. [10] 
were designed in response to NASA’s N+3 initiative to drastically reduce fuel burn, noise, 
and emissions on a 737-800 class airplane with entry-into-service in the 2035 time frame. 
Incorporating a “double bubble” fuselage, the D8 series is a twin aisle design in a 2x4x2 
passenger seating arrangement that takes advantage of a lifting nose reducing the pitching 
moment to trim. The engines, mounted in wing pylons on the 737-800, have been embedded 
in the aft section of the double bubble fuselage. A Pi-tail (7r-tail) arrangement surrounds the 
embedded engines, providing noise shielding while taking advantage of reduced structural 
weight penalties typically associated with a T-tail empennage. The double mounting point 
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of the all-moving horizontal stabilizer to the top of the twin verticals allows for reduced 
structure in both the horizontal and verticals. 

The D8 configuration was sized for a design range of 3,000 nautical miles with a fuel 
reserve equal to 5% of the total trip fuel. To reduce fuel burn, the cruise Mach number 
was reduced from the Boeing 737-800 reference Mach number of 0.78 to Mach 0.72. While 
keeping Mach number constant, the altitude was allowed to vary to minimize fuel burn 
throughout the cruise segment. The climb segment was optimized to minimize fuel burn 
when climbing to altitude, and descent at maximum lift to drag ratio was used. Figure 3.1 
shows a summary of the mission profile used in Ref. [10] and in the multidisciplinary design 
optimization of this research. The mission segments and the different schedules options for 
each are described in Section 5.2.6. 


Cruise 



Figure 3.1: Mission profile used in sizing the D8 concept and the multidisciplinary design 
optimization. 

Transitioned from a Boeing 737-800 design, Ref. [10] describes the incremental aircraft 
changes resulting in the D8.1, a current-technology, double bubble, tube-and-wing concept. 
Modifications included replacing the traditional tube fuselage with the double bubble style 
fuselage, pictured in Fig. 3.2, aft embedded engines, reduced wing sweep from reduced cruise 
Mach number, implementation a doubly supported horizontal stabilizer with the Pi-tail, and 
increased aspect ratio and wing span. The potential benefits of the double bubble fuselage, 
embedded engines, and Pi-tail, were fully described in Ref. [10]. The cruise Mach number 
was decreased from 0.80 (737-800) to 0.72, allowing for reduced wing sweep, increased max- 
imum lift coefficient at takeoff conditions, and reduced structural wing weight. Optimized 
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(a) Fuselage cross-section (b) Fuselage cross-section 

forward of wing box. at wing box. 


Figure 3.2: D8 fuselage cross section forward of the wing box (a) and at the wing box (b) [10]. 

for minimal fuel burn in Massachusetts Institute of Technology’s (MIT’s) conceptual design 
software, TASOPT, the D8.1 three-view is shown in Fig. 3.3 with some geometry parameters 
provided. A more encompassing list of geometric parameters is provided in Table 3.1. Nu- 
merical values with a tilde were not specified in Ref. [10] and had to be inferred, or estimated, 
from drawings provided in the report. 



Figure 3.3: Three view of D8.1 geometry including sectional views [10]. 

With a wing span of 150 feet, the D8.1 configuration benefits were masked by the increase 
in span efficiency from the large span. To better capture the benefits of the double bubble 
configuration, the span was limited to 118 feet, placing it in the same operational class as the 
737-800. This allowed for the capture of the double bubble configuration benefits without 
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Table 3.1: D8.1 geometry parameters from Ref. [10] or as measured from a scaled drawing 
in AutoCAD. 


Parameter 

Symbol 

Value 

Units 

Comment 

Mach number 

M 

0.72 

- 

Start of cruise Mach number 

Wing Area 

S 

1298 

ft 2 

Reference area 

Wing Span 

b 

150 

ft 

Projected wing span 

MAC 

c 

10.6 

ft 

Mean aerodynamic chord 

Aspect Ratio 

AR 

17.3 

- 

- 

Lift-to-Drag 

L/D 

22.1 

- 

Start of climb L/D 

Sweep 

•A-c/4 

—7.2 

deg 

Quarter-chord sweep angle 

Dihedral 

r 

—3.3 

deg 

- 

Taper ratio 

A 

-0.15 

- 

Root chord at fuselage/ wing joint 


the induced drag benefits of a large wing span, which can be added to any conventional 
configuration. The twin engine, span restricted, current technology configuration was called 
the D8.2b. Modifications from the original D8.1 to the D8.2b are described in Table 3.2, as 
optimized in TASOPT for minimum fuel burn. Approximate values were taken from Fig. 3.4. 
Not indicated in Table 3.2 was the change from three rear embedded engines to two engines, 
shown if Fig. 3.4. An Open Vehicle Sketch Pad (OpenVSP) model was created, shown in 
Fig. 3.5, where the engines have been ignored for model simplicity. 

Table 3.2: Parametric geometry changes from D8.1 to D8.2b configuration 


Parameter 

D8.1 

D8.2b 

Wing Area 

1298 ft 2 

1110 ft 2 

Wing Span 

150 ft 

118 ft 

MAC 

10.6 ft 

11.03 ft 

AR 

17.3 

12.5 

L/D 

22.1 

19.45 

Sweep 

—7.2 deg 

—13 deg 

Dihedral 

—3.3 deg 

—5.0 deg 

Taper 

-0.15 

-0.18 


TASOPT predicts a maximum gross takeoff weight of 121,295 pounds with a fuel weight 
of 21,420 pounds. This was the amount of fuel predicted by TASOPT to complete a 3,000 
nm mission with 5% of total fuel reserve. Measurements from Fig. 3.4 placed the fore and 
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Figure 3.4: Three view of the D8.2b geometry, including sectional views. 



Figure 3.5: OpenVSP model of the D8.2b 
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aft limits of the center of gravity to be -12% and 54% mean aerodynamic chord, with the 
MAC specified in Table 3.2. The leading edge of the mean aerodynamic chord was 58.2 feet 
aft of the aircraft nose. A 10% minimum static margin specified the rear CG limit and the 
forward CG limit was determined from a landing/decent condition [10]. It should be noted 
that the elevator deflection angle for the forward CG limit at the landing condition exceeded 
the allowable deflection used in this research and required an adjustment in the forward CG 
location, reduced to 30% MAC, discussed in Section 4.4.2. 

The Flight Optimization System (FLOPS) [14], described in Section 5.2, calculated 
weight estimates that were used in this research instead of inputting the fixed gross weight 
given in Fig 3.4. FLOPS does not require, nor estimate, mass moments of inertia in its 
analysis, but mass moments of inertia were required for the flight dynamic analysis. A 
detailed, high-order analysis for predicting the mass moments of inertia was beyond the 
scope of this research and a simple, conceptual design level, prediction technique was used. 
Radii of gyration, described by Raymer in Ref. [33], were used to estimate the mass moments 
of inertia for the D8.2b geometry. Equation 3.1 gives the inertia calculations in the body 
axis 


h 


h 


b 2 w 


4 g 


-Ri 


L 2 rW 


4 g 


-Ri 



w 


-Ri 


(3.1) 


Depending on the aircraft configuration, the radii of gyration vary requiring the selection of 
a base configuration type to calculate the mass moments of inertia. Reference [33] provided 
a list of aircraft type and the corresponding radii of gyration. A transport aircraft with 
fuselage mounted engines was chosen, and the selected radii are given in Table 3.3. 
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Table 3.3: Radii of gyration used in the mass moments of inertia calculations [33]. 


Radii of Gyration 

Rx 

Ry 

Rz 

Value 

0.24 

0.36 

0.44 


3.2 Cessna 182T Skylane Geometry 


Reference areas and mass properties for the Cessna 182T came from Refs. [15,34], while 
the dimensions and component placements came from the Pilot’s Operating Handbook [16] 
and Specification and Description handbook [35]. Table 3.4 summarizes the reference areas 
and mass properties of the Cessna 182T. Any dimensions not specifically listed in the draw- 
ings of the POH or S&D were measured in AutoCAD 2014 using scaled drawings. Figures 3.7 
and 3.6 were taken from the S&D and POH respectively and scaled, ensuring proper draw- 
ing size and aspect ratio. Using these drawings, the input hies for the aerodynamic analysis 
were created, discussed in Section 5.1. The center of gravity was located at 26.4% from the 
leading edge of the mean aerodynamic chord (c), which was located 25.98 inches from the 
reference datum [16]. The datum — the front face, lower portion of the firewall — was 64.7 
inches measured from the front of the propeller spinner. 

Table 3.4: Cessna 182T reference parameters and mass properties [15]. 


Parameter 

Symbol 

Value 

Unit 

Reference Area 

S 

174 

ft 2 

Mean Aerodynamic Chord 

c 

4.9 

ft 

Wing Span 

b 

36 

ft 

Weight 

W 

2,650 

lb 

x-axis Inertia 

lx 

948 

slug-ft 2 

y-axis Inertia 

ly 

1,346 

slug-ft 2 

z-axis Inertia 

h 

1,967 

slug-ft 2 

xz-axis Inertia 

^xz 

0 

slug-ft 2 


The wing had a span was 36 feet with a wing reference area of 174 square feet, located 
at the top of the fuselage approximately 30 inches above the spinner centerline and 89 inches 
aft of the spinner tip, or 25.98 inches aft of the reference datum. From root to tip, the 
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quarter-chord sweep angle is zero degrees. Broken into two segments, the inboard section is 
of constant chord 64 inches in length, consisting of the NACA 2412 airfoil, with a dihedral 
angle of 2 degrees. The inboard section ends where the strut joins the wing, 46% semispan, 
and the outboard section has a taper ratio of 0.67 and a dihedral angle of 2.5 degrees. The 
ailerons were a constant 8.75 inches in chord starting at the outboard section and extending 
to 97% semispan. A front and top view of the main wing can be seen in Fig. 3.6. 

A conventional tail configuration, the horizontal stabilizer was located 22.2 feet aft and 
1.4 feet above the spinner as extended to the aircraft centerline. Consisting of only one 
trapezoidal segment, the total span was 11.75 feet with a root chord of 4.34 feet, again as 
extended to the aircraft centerline. The taper ratio was 0.62 with no dihedral and a quarter- 
chord sweep angle of 2.25 degrees. A hinge line perpendicular to the x-z plane of the aircraft 
was used for the elevator resulting in a variable chord elevator. At the root, the elevator was 
42.3% local chord and at the tip it was 33.7% local chord. 

Figure 3.7 shows the best view of the vertical stabilizer including the vertical strake 
that extends forward at the base. The vertical stabilizer spans 59.9 inches as measured from 
the top of the fuselage where the vertical the upper surface of the fuselage. At this location, 
the reference root chord (trapezoidal vertical without the strake) was 53.9 inches and the 
tip chord was 27.2 inches. The trapezoidal vertical stabilizer leading edge sweep angle was 
41 degrees and the strake leading edge sweep was 79 degrees. The leading edge of the strake 
joins the upper surface of the fuselage 17 feet behind the front of the spinner. The rudder 
spans the entire vertical stabilizer trailing edge and was nearly constant percent local chord; 
at the tip the rudder was 40% local chord and at the base was 41.7% of the vertical reference 
chord of 53.9 inches. 
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Figure 3.6: Top and front view of Cessna 182T (not to scale) [35]. 
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Fuselage Station (FS) Inches 


Figure 3.7: Side view of Cessna 182T indicating location of mean aerodynamic chord (MAC) 
and datum point (not to scale) [16]. 
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Chapter 4 

System Flight Dynamics and Disturbances 


Described is the stability and control analysis used in the multidisciplinary design opti- 
mization. Figure 4.1 shows the general process for the dynamic analysis with system distur- 
bances. The tan boxes are discussed explicitly with the derivations of the system dynamics 
and the linear quadratic regulator controller presented in Sections 4.1 and 4.2, respectively. 
The continuous turbulence and discrete gust atmospheric disturbances are described in Sec- 
tion 4.3, and the system perturbations, static trim limits, and dynamic response constraints 
are discussed in Section 4.5. Additionally, the flight conditions used for the static trim and 
dynamic analysis are described Section 4.4. Gray boxes with dashed arrows indicate that the 
analysis does not explicitly occur during the dynamic analysis, and the MATLAB® linear 
simulation function was used for all dynamic simulations. 



Figure 4.1: System flight dynamics and disturbances process chart. 
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4.1 Equations of Motion 


The dynamics of the aircraft were modeled by starting with the fully coupled equations 
of motion. These equations were derived and linearized about a steady-state condition 
resulting in a state-space representation of the perturbation equations. State feedback was 
chosen for the controller structure due to the guaranteed closed-loop stability properties 
of the controller, which worked well in a diverse design space explored by an optimizer. 
Although it is rare to have all the states available, the focus of this research was not to 
design a robust controller but rather to incorporate active control into the conceptual design 
space that would give a stable solution. As such, it was undesirable to have the optimizer 
throw out feasible designs due to instability resulting from a poorly designed controller, and 
so a controller design with guaranteed closed-loop stability was used. 

The perturbation equations were derived assuming a symmetric geometry with zero 
steady sideslip. The steady-state thrust terms were solved for explicitly and substituted into 
the state-space form of the perturbation equations. The implicit, state-space form of the 
fully coupled perturbation equations is given by 


Ex = Ax + 


B 


B n 



(4.1) 


where B g is the gust input matrix, and the hat ( ' ) indicates the state and control matrices 
have yet to be premultiplied by the generalized inertial matrix, E. The state and control 
vectors are 


x = 


5 = 


ILg ~ 


uvwpqr(j)6ip 

r 

S e 5 a 5 r 


V* gust Vgust Wgust 


(4.2) 
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with the full definition of the matrices in Eq. 4.1 given in Appendix A. Inverting the gener- 
alized inertial matrix, E, one obtains the standard, explicit, state-space model 


x = E-'A + E- 1 


B 


B n 



Ax + 


B 


B n 



(4.3) 


A linear actuator model was added to the system dynamics for each control surface and 
was modeled by a simple-lag filter transfer function given by [4] 


S(s) _ 1 

u(s) TS + 1 


(4.4) 


where r is the time constant of the filter. The time constant was chosen to be the same 
for all three actuators and was selected as r — 1/(20. 2) s [4], The perturbation equations 
described by Eqs. 4. 1-4. 3 were augmented to include the actuator dynamics while neglecting 
all unsteady terms, thrust terms (except Ct x ), steady-state roll (. R ), and negligibly small 
terms. A vortex lattice aerodynamic analysis tool called Athena Vortex Lattice (AVL), 1 de- 
scribed in Section 5.1, was used for the aerodynamic model, which does not include thrust in 
the analysis, therefore those terms had to be neglected. The velocity dependent thrust term, 
Ct Xu , was included in the perturbation equations as it would be the largest in magnitude for 
the configurations studied in this research. It was approximated using Ref. [32] to be 


Propeller Aircraft: Ct Xu = —3 (C'o + Cw 0 sin#) 
Jet Aircraft: Ct Xu = —2 (Co + CV 0 sin#) 


(4.5) 


The augmented perturbation equations are presented in compact form in Eq. 4.6 with the 
matrix definitions given in Appendix A. Again, the hat used in the Appendix A equa- 
tions indicates that a matrix has not been premultiplied by the inverted generalized inertial 

^ttp : //web .mit . edu/drela/Publi c/web/ avl/ 
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matrix, E aug . 



The actuator input vector, u, is defined as 


u = 


U P Un 


U r 


(4.6) 


(4.7) 


4.2 State Feedback by Linear Quadratic Regulator 

A complete dynamic system of an aircraft incorporating active control is complex, re- 
quiring numerous loop closures to provide adequate closed-loop system response. Classical 
control relics on the iterative selection of gains to achieve the closed-loop system stability, 
but there is no guarantee the gains chosen will be optimum. Modern control theory takes 
advantage of current computing power where numerous linear equations can be solved simul- 
taneously to obtain a set of gains that minimizes a chosen performance index (PI). It is in 
the selection of the performance index that the true engineering of the controller occurs [4] . 

A linear quadratic regulator (LQR) can be used to simultaneously close all the loops in 
a linear, time-invariant (LTI), multi-input multi-output (MIMO) system. With the closure 
of all the loops, the gains are solved for simultaneously, negating the need for successive 
loop closure as required in classical control theory. Extremely versatile, the LQR is capable 
of using performance indices with state and control weighting, time weighting, and deriva- 
tive weighting of the states in both state feedback and output feedback control structures. 
Without any restrictions on the gains, other than closed-loop stability is required, the LQR. 
may choose to zero a gain, thus leaving a loop unclosed. Additionally, compensators may 
be used in the form of filters, integral and derivative controllers. This flexibility makes it an 
excellent tool for finding optimal gains for a controller in an LTI system. 

As a regulator, any non-zero states are driven to zero in such a way that a chosen 
performance index is minimized. By driving all the states to zero, the system is returned 
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to the steady-state condition with the least amount of cost, thus being an optimal solution. 
This is ideal for a stability augmentation system where any deviation from the steady state 
is undesired. 

4.2.1 The Standard Linear Quadratic Regulator Performance Index 

A set of nonlinear equations of motion can be linearized about a condition, a steady-state 
condition for example, and represented in state-space form as given by 

x = Ax + Bu (4.8) 

where both x and u are functions of time. Equation 4.8 is a simplified version of Eq. 4.6. The 
state vector is a vector of perturbations from the steady state condition that the regulator 
drives to zero. With a state-feedback control law 

u = —Kx (4.9) 


the closed-loop system takes on the form 

x — (A — BK) x = A c x (4.10) 

A performance output can be defined as 


z = Hx (4-11) 

where z is a vector of states or a combination of states. An example would be a normal load 
factor as a function of several longitudinal states such as pitch rate, pitch angle, and angle 
of attack. Additionally, for a regulator z can be set to the error of a specific state such as a 
non-zero pitch rate that should be driven to zero. 
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The LQR finds the optimal gains through the minimization of a performance index that 
integrates the values of both the state and control vectors over time. Each state and control 
is weighted in the performance index through the use of weighting matrices, varying the 
impact of each state and control on the performance. The standard performance index for 
the LQR is 

J=- / (x t Qx + u t Ru) dt (4-12) 

2 J o 

with Q > 0,-R > 0. The performance output, z, can be incorporated into the performance 
index by setting Q = H T H. Substituting Eq. 4.9 into Eq. 4.12, the performance index 
becomes 


1 f°° 

J — - ( x t Qx + x t K t RKx ) dt 

° (4-13) 

1 roo v 7 

= - / [a; T (Q + K t RK) x] dt 

2 Jo 


As described by Stevens and Lewis [4] , this dynamic optimization problem can be converted 
into a static optimization problem that is easier to solve. Suppose that a constant, symmetric, 
positive-semidefinite matrix P can be found such that 


( x T Px ) = -x T (Q + K t RK) X 

CLL 


(4.14) 


Substitute the left side of Eq. 4.14 into Eq. 4.13 to get 



Evaluating Eq. 4.15 at the limits of integration gives 

J = ^x T (0)Pa;(0) — ^ lim x^ifjPx^) 

2 2 J — S'-OO 

28 


(4.15) 


(4.16) 



Assuming the closed-loop system is asymptotically stable, x 1 (t)Px(t) will vanish with time 
leaving 

J = i x T (0)Pa;(0) = hr(PX) (4.17) 

where 

X = x(0)a; T (0) (4.18) 

If a positive, semidehnite solution P exists that satisfies Eq. 4.14, a constraint equation can 
be formed by substituting Eq. 4.10 into Eq. 4.14 giving 

-x T (Q + K t RK) x = j f ( x T Px ) 

= x T Px + x T Px (4.19) 

= x T {A t c P + PA c )x 

This constraint equation must be true for all values of x(t) which leaves 

g = A^P+PA C + K t RK + Q = 0 (4.20) 

The time-dependent optimization problem becomes a time-independent optimization prob- 
lem subject to the constraint of Eq. 4.20. 

One can solve a modified problem using the Lagrange multiplier approach which adjoins 
the constraint equation (Eq. 4.20) to the performance index using the Hamiltonian 

J4? = tr(PX) + tr(gS) (4.21) 

where S is a symmetric matrix that needs to be determined, simplifying the constrained 
optimization problem to minimizing Eq. 4.21 without constraints. This is accomplished by 
setting the partial derivatives of the Hamiltonian with respect to each independent variable 
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to zero as follows: 


(4.22) 



the algebraic Riccati equation is found. Because the Kalman gain is independent of S it does 
not depend on the initial state, X. Solving the Riccati equation for the symmetric positive 
semidehnite matrix P, the optimal gain can be found directly using Eq. 4.25. 

4.2.2 Derivation of the Modified Performance Index 

A limitation of the linear quadratic method is the n x n entities that must be chosen in 
the weighting matrix Q where the values may not directly correspond to a performance ob- 
jective due to an observability requirement, initially presented by Kalman [36] and discussed 
by Stevens and Lewis [4], This results in a trial-and-error method of selecting Q where the 
entries are varied until an acceptable transient response is obtained. This method of design 
is highly undesirable. 

Eliminating any restriction of observability in the selection of Q allows for the ele- 
ments to be chosen strictly on desired performance objectives. With specified performance 
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objectives, the structure of the PI and the number of entities to be chosen for the weight- 
ing matrices can be reduced with the closed-loop response dependent on the design of the 
performance index. 

A strength of the linear quadratic method is the selection flexibility of the performance 
index structure where the standard PI, given in Eq. 4.12, can be modified by adding time and 
derivative weighting of the states. A benefit of the time-weighted performance index is that 
it can satisfy the observability requirement, freeing the selection of the weighting matrices. 
The standard PI only lightly penalizes small errors due to a slow pole(s) with small residue, 
resulting in the time to reach a steady state condition being rather large. The time-weighted 
performance index heavily penalizes errors that occur late in the response and, as a result, 
suppresses the effect of a slow pole and lightly damped settling behavior [4], 

The derivative of the state can sometimes be a more accurate representation of the 
workload on a control system and should be weighted in the performance index rather than 
the state itself. For example, the rate of change of a control surface is a more accurate 
representation of required actuator power than the deflection angle itself. 

The standard PI for the linear quadratic regulator was modified to include time weight- 
ing on the states and derivative of the states is given by 

1 f°° 

J=- / (t k x T Px + u t Ru + z T Wz) dt (4.27) 

2 Jo 

where R has the same definition as in the standard performance index, Q has been set equal 
to zero, P is a state weighting matrix for the time weighted terms, and W is a new weighting 
matrix on the states’ time rate of change. Equation 4.27 can be rewritten as 


J=- f (t k e T e + u T Ru + z T W z) dt 

2 Jo 


where e is the error, defined as 


e = z = Hx 


(4.28) 


(4.29) 
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and P = H T H. The solution to Eq. 4.28 can be shown to be 


J= l -tr(P k X ) 

subject to the nested Lyapunov equations given by 


(4.30) 


0 — go = AJPq + PqA c + P 
0 = g\ = A]; Pi + P\A C + P 0 

; (4-si) 

0 = g k - 1 = AjP k _i + P k _\A c + P k -2 
0 = g k = AjP k + P k A c + k\P k _i + + A t c H t WHA c 

with X defined as 

X = x(0)a; T (0) = P'^(0)^ T (0) (P _1 ) T (4.32) 

Derivations of Eqs. 4.30-4.32 are given in Appendix C. Any non-negative integer can be 
chosen for the time-weighting exponent, k, where k = 0 will give the standard performance 
index equations without nested Lyapunov equations. Experience with Ref. [4] indicates that 
k — 2 is a good selection and was used in this research. 

The solution to this minimization problem is dependent upon the initial conditions, x(0), 
unlike the state feedback LQR with the standard performance index. A simple solution for 
eliminating the dependence on the initial state is to average the performance of a set of 
linearly independent initial conditions. This is equivalent to a random variable uniformly 
distributed on the surface of a unit sphere of dimension equal to the length of the state 
vector [37]. In essence, instead of minimizing the performance index in Eq. 4.30, the expected 
value of J, E{J}, is minimized such that 

E{J} = ^E{x T (0)P k x(0)} = l -tr(P k X ) (4.33) 
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where the symmetric n x n matrix 


X = A{a;(0)a; T (0)} (4.34) 

is the initial autocorrelation of the state. For the regulator problem, it is practical to set 
X = I since it is desired to drive arbitrary nonzero states to zero [4]. 

4.2.3 Selection of the Weighting Matrices 

An emphasis was placed on reducing the number of entities that must be chosen in the 
weighting matrices. As a regulator with state feedback, any nonzero state was considered 
an error. The performance output, z, was selected as 

r i t 

z = e = Hx =ul3apqr(j)0'ip8 e 5 a 8 r (4.35) 

with the performance output matrix H defined as 

H =™ dia s{ jt a i i 1 1 1 1 1 1 1 1 l} (4.36) 

The performance output was chosen to give one degree of error equal weighting as one foot 
per second in the performance index. This was chosen to maintain the generality of the 
controller while reducing the number of entities to be chosen in the performance index. 

The derivative- weighting matrix W was chosen as 

W = diag jo 0000000011 lj (4.37) 

to place a penalty on high control surface deflection rates. 
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The controls-weighting matrix R was selected as 


1 0 0 


R = 0.1 


0 1 0 


0 0 1 


(4.38) 


where smaller values of R give greater authority to the control inputs. However, values less 
than those given in Eq. 4.38 produce transient responses negligibly different than results 
using Eq. 4.38. 


4.3 Atmospheric Disturbances 

To stress the control system, atmospheric disturbances were modeled and used to test 
the response of the system to deviations from the steady-state condition. Including these 
disturbances adds validity to the model and eliminates unrealistic configurations. Two atmo- 
spheric disturbance models were used as suggested in military standard MIL-STD-1797A [9]: 
a continuous turbulence model in the frequency domain and a discrete gust model in the 
time domain. 


4.3.1 Continuous Turbulence 

The von Karman continuous turbulence model was used as specified in MIL-STD-1797A. 
The one-sided spectra are of the form [29,30,38] 




2"l 5/6 


[l + (1.339L„fi) 2 ] 

,2 L, 1 + | (2.678L„Si) 2 

^v g \^ L ) u v 11/6 

* [1 + (2.678L„fi) ] ' 

* 1 + § (2.678U,!2) 2 

. otII/6 

71 [1 + (2.678 L W Q) 2 ' 


(4.39) 
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where L n is the length scale and a n is the root-mean-square intensity of the continuous 
turbulence, n = u,v,w. The spatial frequency, used in the von Karman form of the 
spectra is related to the temporal frequency, oj, by = oj/U. The spectra are converted to 
functions of u using the simple relationship 

$n 9 M = n = u,v,w (4.40) 

The mean-squared response of each state can be found be integrating the magnitude-squared 
system transfer function, |G(hn)| 2 , times the power spectral density of the continuous tur- 
bulence response at all frequencies, shown in Eq. 4.41. 

/»oo 

a l n = \G(iuj)\ 2 $ ng (uj)duj, n — u,v,w (4.41) 

Jo 

Medium to High Altitude Clear Air Turbulence 

Isotropic turbulence is assumed for clear air turbulence at altitudes of 2,000 feet and 
above. The turbulence length scales and mean-square intensities are defined as [29, 30, 38] 


2 2 2 
<7Z = at = err, 


L„ = 2 L„ = 2 L„ 


(4.42) 


The length scale for medium to high altitude turbulence is L u = 2, 500 ft for the von Karman 
form and L u = 1,750 feet for the Dryden form, which is used in the discrete gust model 
described in Section 4.3.2. 

The RMS intensities to be used in Eq. 4.39 were specified in Table 4 of Section 3. 2. 1.5.1 
of SAE-AS94900 [39]. The table gives the vertical RMS intensities for numerous altitudes 
with varying probability of exceedance, and the longitudinal and lateral intensities can be 
found by Eq. 4.42. Light, moderate, and severe levels, with probabilities of exceedance of 
10 -2 , 10~ 3 , and 10~ 5 , are given in Table 4.1 and follow the turbulence severity levels shown 
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in Fig. 4.2. Altitudes used between data points were linearly interpolated to allow for any 
desired altitude between 2,000 and 85,000 feet. 


Table 4.1: RMS gust intensities, fps. [39]. 


Altitude (ft) 

Light 

Moderate 

Severe 

500 

6.6 

8.6 

15.6 

1,750 

6.9 

9.6 

17.6 

3,750 

7.4 

10.6 

23.0 

7,500 

6.7 

10.1 

23.6 

15,000 

4.6 

8.0 

22.1 

25,000 

2.7 

6.6 

20.0 

35,000 

0.4 

5.0 

16.0 

45,000 

0 

4.2 

15.1 

55,000 

0 

2.7 

12.1 

65,000 

0 

0 

7.9 

75,000 

0 

0 

6.2 

85,000 

0 

0 

5.1 
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Figure 4.2: Turbulence severity and exeedance probability [9]. 

Low Altitude Clear Air Turbulence 

Isotropic turbulence can no longer be assumed at altitudes below 2,000 feet, invalidating 
Eq. 4.42. In this altitude region, the von Karman and Dryden models use the same length 
and RMS intensity, depending on altitude and probability of exceedance. This altitude 
region is broken into two regions: less than or equal to 1,000 feet, and between 1,000 feet 
and 2,000 feet. Nonlinear relationships between the turbulence length and intensity scales 
are used in the lower region altitudes. For the intensity scales, the longitudinal and lateral 
intensities are related to the vertical intensity scale by 


O U O V 


(0.177 + 0.000823/r)°' 4 


(4.43) 
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where h is altitude and a w is specified by 


<7 W — 0. Ilt20 


(4.44) 


The term w 2 o is the mean wind speed measured 20 feet above the ground, a function the 
intensity probability of exeedance. Figure 4.3 gives the means wind speed 20 feet above 
ground level for varying probabilities of exceedance. The intensities used for light, moderate, 
and severe were 23, 30, and 45 knots respectively; these correlate to probabilities of exeedance 
of 10" 2 , 10~ 3 , and 10~ 5 . 

The turbulence length scales are given by 


2 L w = h 
L u 


h 

(0.177 + 0.000823h) L2 


and, as with the intensities, h is altitude. 


(4.45) 
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Power-series Approximation of Linear-gradient Turbulence 


In addition to uniform immersion of the center of gravity in turbulence, captured by u g , 
v g , and w g , linear gradients of the turbulence velocities can be modeled [8,9,29,30,40-42], 
Linear gradients of the turbulence held are, in effect, equivalent to aircraft angular velocities 
and can be modeled as such. Multiple notations are used for the linear gradients (Ref. [43]) 
but the forms given in MIL-F-8785C were chosen to be presented here as the associated 
spectrums were also given. The equivalent angular velocities from variations in the gust held 
across the length and span of the aircraft are defined as [8] 


Pa 


-a = q g 
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_dw l 

dy 

dw a 

dx 

_dv l 

dx 


The spectrums associated with Eq. 4.46 are 


(4.46) 
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(4.47) 


where $ Ws and were dehned in Eq. 4.39. To remain consistent with the turbulence length 
scale definitions given in Eqs. 4.42 and 4.45, was modified from MIL-F-8785C by adding 
a factor of two to the L w term. 

Inclusion of this analysis was considered for this research, but it was decided not to 
implement it into the current methodology. Reviewing Refs. [8,9,29,30,39,44], the rotations 
are added as another input with three new spectrums, given in Eq. 4.47, and therefore 
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produce a separate set of mean-squared responses, responses that do not have an associated 
requirement in SAE-AS94900. SAE-AS94900 only gives requirements for the RMS gust 
response to u, v, and w gusts at the center of gravity. Another limitation of the linear- 
gradient approximation is the limited range of validity of the approximation; only at low 
frequencies does the approximation apply, wavelengths on the order of eight times the total 
vehicle length and greater. In the case of the discrete gust, discussed in Section 4.3.2, 
the gust is tuned to the natural frequencies of the closed-loop system which change with 
every iteration. This means the approximation could be valid in some cases while not 
in others, making the analysis inconsistent across all configurations. Finally, structural 
modes analysis was not included in this research and therefore the current proposed analysis 
provides a sufficient dynamic constraint for incorporating flight dynamics and control into 
the conceptual design process. 

4.3.2 Discrete Gust 

The discrete gust model has the “1-cosine” profile defined by Eq. 4.48 and illustrated 
in Fig. 4.4 


v = 0, 





x < 0 

0 < x < d m 


v = 0, 


X > d n 


(4.48) 


where V m is the magnitude of the gust and d m is one-half the total gust length. The gust 
length was calculated to provide maximum system excitation by tuning the gust to the least 
damped, one longitudinal and one lateral, system natural frequencies [39]. The magnitude 
of the gust, V m , was calculated using Fig. 4.5 [9]; P(0.01) corresponds to light turbulence 
intensities and P(0.001) corresponds to moderate turbulence intensities. The Dryden length 
scales were used with Fig. 4.5 as directed by MIL-STD-1797A for the discrete gust model, and 
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Figure 4.4: “1-cos” discrete gust profile [30]. 


SAE-A94900 specifies the discrete gust peak amplitude of 60 fps for severe turbulence. The 
moderate turbulence intensity curve in Fig. 4.5 was added in addition to the light turbulence 
intensity curve given in MIL-STD-1797A [9]. The data used to create the second curve was 
presented by Roskam in Ref. [30] and is summarized in Fig. 4.6. The analysis is capable of 
using either the calculated gust magnitude taken from Fig. 4.5, or the severe turbulence case 
of 60 fps. 



■o 

<u 

N 

W 


£ 

e 

> 


© 

3 


<D 

■o 


C 

O) 

rc 


w 

3 

0 



Normalized Discrete Gust, 7 

L 2L 2L 

U V w 


Figure 4.5: Normalized discrete gust for determining gust magnitude [8,9]. 
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Figure 4.6: Probability of equaling or exceeding a given gust magnitude [30]. 

4.4 Flight Conditions for Dynamic Response Analysis 

A single flight condition cannot be used to fully analyze the capability of a flight control 
system. Numerous flight conditions must be used to ensure the flight control system is able to 
provide acceptable dynamic response throughout the flight envelope for a given configuration. 
Both nominal and degraded flight control system capabilities directly affect the sizing of 
the stabilizer and control surfaces. In this research, only nominal flight control system 
performance was explored as bounds for the conceptual design space were desired without 
the necessity of robust control system design. Once a configuration is chosen for preliminary 
or detailed analysis, a robust controller, including failure modes, must be designed to ensure 
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adequate handling qualities. This detailed controller design was beyond the scope of this 
research. 

4.4.1 Cessna 182T Flight Condition for Tool Development and Exercise 

The Cessna 182 is a very common four place, single-engine aircraft. A model of this 
airplane was used to verify that the various analysis components of the flight dynamics 
analysis tool were working as expected before proceeding with the full system analysis. As 
such, only a single flight condition was used. Experimental data were available for the Cessna 
182T at the cruise condition from Refs. [15,34], At 5,000 feet above sea level, the cruise 
velocity was 220.1 fps and the Mach number was 0.201. This gives a dynamic pressure of 
49.6 lbs/ ft 2 and an air density of 0.002048 slug/ft 3 . This, in combination with the Cessna 
geometry data from Section 3.2 was used as inputs for the stability augmentation system. 

4.4.2 Flight Conditions for Evaluating D8.2b Closed-loop Performance 

Adequate sizing of the control effectors is essential to the controllability and safety of 
any aircraft configuration despite their detrimental effect on system performance. Identifying 
the design constraining flight conditions for sizing the control effectors in conceptual design 
is key to linking the conceptual phase and flight test phase of a design. It is undesirable to 
evaluate the design stability and control characteristics at every point in the flight envelope, 
whereas exploring a few specific, worst-case conditions that constrain the control effectors is 
desired. Chudoba and Cook strived to identify in Ref. [45] a set of design- constraining flight 
conditions (DCFC) for sizing control effectors in the conceptual design phase. 

Using the work of Chudoba and Cook in Ref. [45] as a guide, five flight conditions, 
with nominal flight control system performance, were chosen to evaluate the adequate sizing 
of the control effectors and closed-loop dynamic response to system disturbances. Three 
steady-state flight conditions were used to evaluate the control effectors’ ability to properly 
trim the aircraft under longitudinal and lateral/directional loads. At takeoff, a nose up trim 
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condition at rotation speed was used to evaluate the control power of the elevator with the 
center of gravity in the most forward position. At the most forward CG position, the greatest 
amount of pitch trim control power is required, constraining the size of the elevator. The 
methodology used in this research establishes a rear-most allowable center of gravity position 
but may not capture the pitch-trim demands of a takeoff condition. To model this, a center 
of gravity range of 30% mean aerodynamic chord was used to establish the forward-most 
position, based on the CG envelope of the Boeing 737-800 shown in Fig. 4.7. 



Figure 4.7: Center of gravity envelope for Boeing 737-800 aircraft [46]. 


For the second static case, a one engine inoperative (OEI) condition was chosen as a 


test of the lateral/direction control power of the control system. Lateral/directional control 


effectors must be capable of trimming the adverse moments, specifically a yawing moment, 
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created during an engine failure occurring after the takeoff decision speed. In a conventional 
configuration empennage, i.e. the rudder is the main effector for producing yawing moments, 
the OEI takeoff condition can be an active constraint on the sizing of directional control 
effectors. 

The last static case was a low-speed maneuver condition with the center of gravity 
in the aft position. In the clean configuration, a loading condition of 3.75g’s was used to 
calculate a constant pitch rate required of the control system. This was designed to test the 
control system under a heavily loaded flight condition at the clean configuration maximum 
lift coefficient. 

Two dynamic flight conditions, both at the aft-most center of gravity position, were 
selected to evaluate the dynamic response performance of the system to disturbances. A lg 
stall flight condition tests the flight control system in a low dynamic pressure state where 
the disturbance will place a heavy workload on the stability augmentation system, causing 
the dynamic constraints to become active in the design space. The second dynamic flight 
condition evaluated was the cruise condition where the dynamic response was evaluated in 
the mission segment with greatest influence on performance. The altitudes for the stall and 
maneuver conditions were chosen to be 2,000 feet due to the adverse atmospheric condi- 
tions at the lower altitude. For cruise, a target cruise altitude of the D8.2b configuration, 
calculated by the Flight Optimization System, described in Section 5.2, was used. The num- 
ber of dynamic flight conditions was limited to reduce the overall analysis time. Table 4.2 
summarizes the chosen flight conditions. 
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Table 4.2: Flight conditions for D8 system evaluation. 


Analysis 

Type 

Description 

Configuration 

CG 

Position 

Atmospheric 

Disturbances 

Trim 

Nose up with rotation at main 
gear ground contact point (SL) 

Takeoff 

Forward 

No 

Trim 

Single engine out at 35 ft 

Takeoff 

Aft 

No 

Trim 

Maneuver speed at 2,000 ft 

Clean 

Aft 

No 

Dynamic 

lg stall at 2,000 ft 

Clean 

Aft 

Yes 

Dynamic 

Cruise at 40,000 ft 

Clean 

Aft 

Yes 


4.5 Static Trim and Dynamic Response Constraints 


Reducing horizontal and vertical stabilizer areas results in degraded flying qualities of 
any configuration, and the damping of dynamic modes is directly affected by the stabilizer 
areas. Control power is decreased with decreasing empennage size due to reduced control 
effector area, thus requiring flying quality enhancements to be provided through an active 
control stability augmentation system. The control system workload increases with decreas- 
ing empennage area to a point where the control system cannot provide adequate flying 
qualities due to limitations in the control system, such as maximum deflection angles and re- 
duced control effectiveness from viscous effects. To capture these control system constraints, 
maximum control surface deflection angles were chosen and dynamic response requirements 
were specified. Failure to meet these constraints indicates the system was unable to ade- 
quately stabilize the configuration to maintain adequate handling qualities. 


4.5.1 Aileron and Rudder Deflection Limits 

For the lateral/directional control surfaces, ailerons and rudders, a twenty degree max- 
imum control surface deflection angle was chosen as an upper bound for the control system. 
Twenty degrees was chosen to prevent the control system from becoming saturated, removing 
any control authority for commanded controls from the pilot, or autonomous control system, 
while avoiding the nonlinear effects of very large control surface deflections. The total control 
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surface deflection angle was calculated as the sum of the steady-state deflection angle and 
any perturbation deflection angle. This maximum deflection angle was checked during the 
discrete gust simulations, discussed previously, and the steady-state perturbations discussed 
later in this section, to ensure any constraint violation was accounted for, and penalized. 

4.5.2 Cessna 182T Elevator Limit 

A simplified elevator constraint was applied to the Cessna 182T flight condition de- 
scribed in Section 4.4.1. Similar to the aileron and rudder constraints, the elevator was 
limited to a maximum deflection angle of 20 degrees. Again, this was chosen to provide 
additional control authority to the pilot or autonomous system while avoiding the nonlinear 
aerodynamics of flow separation on the control surface. 

4.5.3 D8.2b Elevator Constraint Derivation 

A more comprehensive elevator constraint was appropriate for the D8.2b flight con- 
ditions described in Section 4.4.2. The all-moving horizontal stabilizer has aerodynamic 
properties of a wing rather than the properties of a control surface deflecting aft of a lift- 
ing surface. A deflection angle of 20 degrees, used for the Cessna 182T elevator constraint, 
would be in or beyond the stall region of the D8.2b horizontal stabilizer. Large downwash 
angles would also be present at high lift conditions, reducing the effective angle of attack 
experienced by the horizontal stabilizer. As a result, a geometric deflection angle was not 
appropriate for the application of an D8.2b elevator constraint. 

The horizontal stabilizer airfoil was extracted from the OpenVSP D8.2b horizontal tail 
geometry described in Section 3.1, and shown in Fig. 4.8. The 12% thick airfoil has negative 
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camber of -0.8% providing a download at zero geometric angle of attack. As this custom air- 
foil was used for the horizontal stabilizer, the aerodynamic properties could not be referenced 
from any sources. Therefore, the airfoil aerodynamic properties had to be predicted using an 
analytic method. XFOfL, 2 a subsonic airfoil development system, was used in the prediction 
of the airfoil lift properties, including maximum lift [47-49] . This airfoil design tool is capa- 
ble of viscous analysis with free boundary layer transitions, and maximum lift predictions. 
An angle of attack sweep was performed in XFOfL to determine the airfoil maximum and 
minimum lift coefficients at Mach 0.26 and a Reynolds number of 10.3 million. Mach 0.26 
corresponds to the velocity of the baseline D8.2b configuration at the end of the takeoff roll, 
and beginning of rotation, at sea level, standard atmospheric conditions. The horizontal tail 
mean aerodynamic chord (5.5 feet) was used as the reference length in the Reynolds number 
calculation. Figure 4.9 shows the lift coefficient versus angle of attack as analyzed in XFOfL. 
From Fig. 4.9, the maximum and minimum lift coefficients were determined to be +1.37 and 
-1.50 respectively. 



Angle of Attack, deg 


Figure 4.9: Lift coefficient versus angle of attack for the horizontal tail airfoil at Mach 0.26 
and Reynolds number of 10.3 million. 

The horizontal tail was used for longitudinal pitch control and should not be operated 
at the airfoil stall condition, otherwise a loss of control authority would occur resulting in 
2 http : //web .mit . edu/drela/Public/web/xf oil/ 
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potential loss of control. Therefore, 75% of the section lift coefficient limits were used in 
determining the maximum effective angle of attack constraints for the horizontal stabilizer, 
forcing the stabilizer to remain in the linear range of the lift curve slope. This corresponded 
to section lift coefficient limits of +1.03 and -1.13. 

The baseline horizontal stabilizer, described in Sections 3.1 and 5.1.2, was modeled in 
isolation using Athena Vortex Lattice (AVL), a vortex lattice solver described in Section 5.1. 
The vortex lattice code was used to capture the finite span effects from the tip vortices, and 
to determine the angles of attack where the section lift coefficient equaled either the +1.03 
or -1.13 limit. Since the horizontal stabilizer was an all- moving control surface, the angle of 
attack was varied in AVL, both positive and negative, until the section lift coefficient limit 
was reached. Figure 4.10 shows the horizontal stabilizer modeled in AVL with the section 
lift coefficient equal to +1.03. This occurred at 12.2 degrees angle of attack. 

AVL provides two outputs for the loads on a single lifting surface. The first is to nor- 
malize by the input reference area, typically the wing reference area. The second involves 
calculating the area of each lifting surface and normalizing the load produced by each lifting 
surface by their respective areas. Normalizing the horizontal stabilizer load by its calcu- 
lated area allowed a connection to be made between the full configuration model and the 
horizontal tail modeled in isolation, even if the horizontal tail area changed. As a result 
of the configuration changes to the horizontal tail area during the design optimization, it 
was necessary to use the surface lift coefficient to maintain a connection results similar to 
those of Fig. 4.10. Using the horizontal stabilizer surface lift coefficient made it possible to 
determine if the section lift coefficient was exceeding the allowable limits in the presence of 
downwash/upwash. At the section lift coefficient limits of +1.03 and -1.13, the horizontal 
stabilizer surface lift coefficients of 0.9553 and -1.054 corresponded to angles of attack of 12.2 
and -11.4 degrees. At zero angle of attack, the lift curve slope was calculated in AVL to be 
4.943 per radian. 
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(a) Horizontal stabilizer in AVL. 



Figure 4.10: Simplified D8.2b horizontal stabilizer as modeled in AVL with the associated 
spanwise lift and load distributions. Analysis was run at 12.2 degrees which corresponds to 
maximum section lift coefficient. The dashed line is the section lift coefficient and the solid 
line is the load distribution. 
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As the horizontal tail area varies in the optimization, the aspect ratio and sweep remain 
constant meaning that the lift curve slope was also constant for a specified Mach number. 
With the surface lift coefficient correlation and lift curve slope, it was possible to calculate the 
effective angle of attack and determine if the section lift coefficient limit had been exceeded. 
The effective angle of attack was calculated by 


a UT efJ - 


C L 


surf 


Cr 


+ a HT L= o 


(4.49) 


where ctHT L=0 was the zero lift angle of attack of the horizontal stabilizer, calculated to be 
-1.0 degrees. Using the lift curve slope, surface lift coefficient, and zero lift angle of attack, 
the effective angle of attack could be inversely determined and verified in AVL by setting 
the surface lift coefficient. 

By solving for the effective angle of attack using Eq. 4.49 from the surface lift coefficient, 
a constant constraint could be established for the horizontal stabilizer that would be sensitive 
to the downwash effects at high loading conditions. The horizontal stabilizer effective angle 
of attack constraint was set to —11.4 < QtHT eff < 12.2 degrees. This constraint was used in 
all flight conditions described in Section 4.4.2 for the D8.2b configuration. In the dynamic 
flight conditions, the elevator deflection angle was added to the steady-state effective angle 
of attack to ensure the maximum deflection did not violate the constraint. 

All the flight conditions, with the exception of cruise, were at similar Mach and Reynolds 
numbers, which means the lift curve slope and airfoil maximum lift coefficients were virtu- 
ally constant. In the cruise case, the Mach number was significantly higher at Mach 0.72, 
but the Reynolds number was essentially constant due to the high altitude operation. The 
increase in Mach number reduces the airfoil maximum lift coefficient and increases the lift 
curve slope, which will affect the constraints. However, a constant constraint was required 
for the optimization and therefore the same horizontal tail constraint used in all flight condi- 
tions, resulting in optimistic horizontal stabilizer constraints in the cruise condition. It was 
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predicted, and supported by results shown in Chapter 7, that the cruise condition was not 
constraining on the elevator and was therefore accepted as a constraint. 

4.5.4 Dynamic Response Constraints 

Military specification SAE-AS94900 [39] specified performance requirements on the tran- 
sient response of the attitude Euler angles, 0, 9 , and 0, that allowed for quantitative eval- 
uation of the flying qualities without having to identify specific dynamic modes; this was 
advantageous because the traditional dynamics modes may not exist due to the complex dy- 
namics provided by unconventional configurations, and configurations with an active control 
system. Assessing the transient response allowed for the flying qualities to be evaluated for 
any geometric configuration with any active control system. 

Specified in SAE-AS94900, with active control, the root- mean-square deviations in pitch 
attitude angle, 9 , must be less than or equal to five degrees in a continuous, one-dimensional 
turbulence held. In response to a five degree pitch perturbation, the control system must be 
capable of returning the pitch attitude to within plus or minus 0.5 degrees of the steady-state 
condition within five seconds for aircraft in classes I— III, defined in MIL-STD-1797A [9]. Sim- 
ilar to the pitch attitude, the roll attitude angle, 0, root-mean-square deviation in continuous 
turbulence must be less than ten degrees, and reach a static accuracy of plus or minus one 
degree within five seconds from a five degree roll perturbation. In continuous turbulence, 
the heading angle, 0, must have a root-mean-square heading deviation of less than or equal 
to five degrees [39]. 

While climbing at a maximum rate of 2000 feet per minute, the control system must be 
capable of leveling off and achieving a static airspeed accuracy of plus or minus 10 knots or 
2% of the reference airspeed, whichever is greater. This accuracy must be achieved within 30 
seconds of engaging the airspeed hold. Any residual oscillations within the static accuracy 
margin must have a period of oscillation greater than 20 seconds. This requirement was 
modeled as a small perturbation in the pitch angle, 6 . This dynamic constraint was modeled 
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as a perturbation in pitch where the climb rate was used to calculate the flight path angle, 7 . 
The steady-state angle of attack was then added to the flight path angle, shown in Fig. 4.11, 
to give the total perturbation that the control system must return to within the steady-state 
tolerance in 30 seconds. The dynamic constraints are summarized in Table 4.3. 



Figure 4.11: Calculation of the pitch Euler angle used in the airspeed hold perturbation. 


Table 4.3: Summary of the static trim and dynamic response performance constraints for 
the D8.2b configuration. 


Description 

System Disturbance 

Constraint 

max aileron/rudder deflection 

all 

-20 < 5 < 20 deg 

max elevator deflection 

all 

-11.4 < a HTeff < 12.2 deg 

5 deg pitch 

perturbation 

± 0.5 deg in < 5 s 

5 deg roll 

perturbation 

± 1.0 deg in < 5 s 

airspeed hold 

perturbation 

± 10 kts or 2% < 30 s 

pitch deviation 

co nt. turbulence 

o-e < 5 deg 

roll deviation 

co nt. turbulence 

cr $ < 10 deg 

heading deviation 

cont. turbulence 

cty < 5 deg 
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Chapter 5 
Analysis 

5.1 Aerodynamic Analysis with Athena Vortex Lattice 

The aerodynamic modeling method used in this research was a vortex lattice called 
Athena Vortex Lattice (AVL). Developed and written at MIT, originally by Harold Youngren 
in 1988 and continuously updated by Youngren and Mark Drela since, the methodology was 
based on the original works of Lamar [50-53] , Lan [54] , and Miranda [55] at NASA Langley 
Research Center [56]. A low-order aerodynamics tool, the vortex lattice method was chosen 
due to its low computational cost while providing estimates of the aerodynamic forces and 
moments, capturing the coupled effects between surfaces and limited span. MIT’s AVL 
was chosen due to the code’s ability to run autonomously by reading in a script, its ability 
to simultaneously solve for both the longitudinal and lateral/direction stability derivatives, 
and the open source availability of the code. 1 Additionally, the input and output hies are 
formatted in a way that was easy couple with the other processes in this research. 

Typically when using a vortex lattice tool, only the lifting surfaces are modeled. This 
has been successful in capturing the lift and drag forces on a simple geometry. However, in 
a longitudinal and lateral/directional stability analysis, this representation of the geometry 
was not sufficiently accurate in predicting the aircraft neutral point, a key prediction in the 
implementation of the methodology used in this research. Fuselage effects and propulsion ef- 
fects, destabilizing depending on propulsion type and location, are not captured by modeling 
the lifting surfaces alone and will erroneously predict the neutral point location. Due to the 
low-order nature of the vortex lattice method, special modeling practices must be employed 
to simulate the effects off the propulsion system and the fuselage on aircraft stability. 

ffittp : //web .mit . edu/drela/Publi c/web/ avl/ 
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The fuselage can be simulated as a zero camber lifting surface the shape of the outer 
mold line as viewed from the top and side of the aircraft. Previous work by the author in 
Ref. [57] showed that, for wing-body configurations having either a high- or mid-wing, a 
flat-panel surface approximating the shape of the fuselage, as cut by the horizontal plane 
shown by Fig. 5.1, produced results that matched experimental data well. The original 
geometry from Ref. [58] was a rectangular wing with a symmetric airfoil and a body of 
revolution about the x-axis. The vertical placement of the wing was dependent upon the 
configuration. As can be seen in Fig. 5.1, the body of revolution was modeled as a 2-D 
surface with the shape matching the outer-mold-line as the x-y plane cuts through. The 
high-wing modeling methodology was employed for the Cessna 182T and the mid-wing for 
the initial D8 modeling. 




Figure 5.1: Wing-body configurations modeled in AVL. 


Figures 5.2 and 5.3 show the comparison of the calculated lift and moment curves 
versus angle of attack to the high- wing configuration experimental data. As indicated in the 
figures, AVL predicted the lift curve slope reasonably well, only a small under-prediction 
of the experimental result, while the moment-curve slope was accurately predicted with a 
small offset in the absolute value. In a stability analysis, the slope of the curve is of greater 
importance. 

The mid-wing configuration was modeled differently than the high-wing configuration. 
If the wing was brought all the way to the centerline it would interfere with the fuselage 
lifting surface, giving poor results. To address this issue, the wings were truncated at the 
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Figure 5.2: Comparison of high-wing configuration lift coefficient with experimental data. 

side of the fuselage and the COMPONENT (legacy INDEX) command was used within AVL 
to indicate that the fuselage and wing were joined, properly distributing the vorticity across 
the geometry. This is illustrated in Fig 5.1(b). 

When reexamining the previous work of Ref. [57], a modeling error was discovered for 
the tapered lifting surface that produced erroneous results. A slight overlap existed between 
the wing and the tapered lifting surface that was corrected with the new results presented 
here. In Fig. 5.4, the original lift curve slops is seen to be much larger than the experimental 
results. With the correction of the geometry error, the lift coefficient AVL results now agree 
very well with the experimental results. Surprisingly, the erroneous geometry coincidentally 
happened to produce an accurate estimate of the pitching moment coefficient. With the 
corrected geometry, there was a slight under prediction of pitching moment coefficient slope. 


•••*•• High-wing Exp. Data 
— - AVL Tapered Lifting Surface 
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Figure 5.3: Comparison of high-wing configuration pitching moment coefficient with exper- 
imental data. 
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Figure 5.4: Comparison of mid-wing configuration lift coefficient with experimental data. 
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Figure 5.5: Comparison of mid-wing configuration pitching moment coefficient with experi- 
mental data. 
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5.1.1 Cessna 182T Aerodynamic Modeling 


With the modeling strategies just discussed, a full Cessna 182T AVL geometry was 
created. The complex Cessna 182T geometry described in Section 3.2 was simplified down 
to the lifting surfaces and a cruciform style fuselage. Following the fuselage outer mold line, 
the vertical cross section was at the aircraft plane of symmetry, and the horizontal cross 
section was located in a plane offset three inches below the x-y plane, going through the tip 
of the propeller spinner. The wing struts were not modeled. A rectangular lifting surface 
was placed in front of the geometry at the spinner to model the destabilizing effects of the 
spinning propeller. The horizontal and vertical cuts of the fuselage were modeled with the 
NOWAKE feature of AVL turned on. For those surfaces, no trailing edge wake was shed, 
eliminating any interference of coplanar components downstream, such as the propeller and 
the horizontal fuselage. The NOWAKE feature captured the moments produced by the 
surfaces but no lift was produced. This is a reasonable approximation given that a fuselage 
is a very poor lifting device but it produces a destabilizing effect, both modeled by a fuselage 
without trailing edge vortices. A detailed description of the propulsion effects modeling is 
described in Section 6.3.1. 



Figure 5.6: Cessna 182T geometry as modeled in AVL. 

The Cessna 182T horizontal and vertical stabilizer areas were modified with a fixed 
center of gravity in a simple, two-variable trade study where the flight dynamics and system 
disturbances were tested. The results from this testing are given in Section 6.4. 

In this initial testing, the only geometrical parameters varied were the horizontal and 
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Figure 5.7: Cessna 182T geometry with modified horizontal tail as percentage of baseline 
volume coefficient. 


vertical tail volume coefficients. The vertical tail strake, shown in Fig. 5.6, was neglected to 
maintain the trapezoidal geometry parameters. To reduce the number of variables, the hor- 
izontal and vertical tail moment arms, aspect ratios, quarter-chord sweeps, dihedral angles, 
taper ratios, and airfoils were unchanged from the baseline model. Regardless of the change 
in volume coefficient, the root quarter-chord coordinate locations were fixed for both the 
horizontal and vertical stabilizers. Changing the volume coefficient with these restrictions 
effectively resulted only in a change of reference area and span. 

The horizontal tail volume coefficient is defined as 


C„ = (5.1) 

Traditionally, l u is defined as the length between the wing aerodynamic center and the 
stabilizer aerodynamic center. For simplicity due to the current fixed CG position in this 
analysis, the length Ih used in Eq. 5.1 was measured from the CG position to the root quarter- 
chord point on the stabilizer. Of note, the volume coefficient for the vertical stabilizer is 
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the same equation as Eq. 5.1 with the subscript H replaced by subscript V and mean 
aerodynamic chord, c, replaced by the span, b. 

The tail volume coefficients were varied from the Cessna 182 model values of = 0.7 
and Cy = 0.28 to 25% of the baseline over ten equally spaced steps. This resulted in 100 
geometries evaluated, each with differing stabilizer areas and spans. As the tool was still 
in development, computation time was sacrificed to provide extra data, ensuring it was 
performing as expected. 

Figure 5.7 shows a drawing from the top view of the AVL geometry with varying hor- 
izontal tail volume coefficients. The coefficients shown are 0.70, 0.52, 0.35, and 0.17 which 
correspond to 100%, 75%, 50%, and 25% of the baseline. The centerline of the fuselage is 
indicated by the solid line drawn through the symmetric plane. 

Any experienced aircraft designer would notice the small horizontal tail area for the 25% 
baseline volume coefficient and doubt its ability to achieve adequate handling qualities, much 
less trim at takeoff conditions. As discussed later in Section 7, all cases passed the dynamic 
performance checks for adequate handling qualities. However, focusing only on the cruise 
flight condition was inadequate for accurately sizing the stabilizers, and more constraining 
flight conditions must be included in the analysis. 

To capture any drag benefit from reducing the stabilizer area, and including active 
control in the design process, the reduction in parasitic drag must be calculated from the 
reduced wetted area. A zero-lift drag coefficient in the cruise condition was given by Napoli- 
tano as Cp 0 = 0.027. This is the total parasitic drag for the Cessna 182 Skylane which was 
modeled using AVL. Reducing the horizontal and vertical stabilizer areas will reduce the 
total parasitic skin-friction drag. To model this, the contribution of the baseline empennage 
to the total parasitic drag was estimated using Raymer’s component buildup method where 
each component is denoted by the subscript c [1]. The subsonic parasitic drag coefficient can 
be approximated by 


(Cdo) 


subsonic 


EC fc FF c Q c S wetc 

S 


(5.2) 
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where Cf is the skin- friction drag coefficient, FF is a form factor that estimates the pressure 
drag due to viscous separation, Q is an interference factor, and S wet is the estimated wetted 
area of the component. The miscellaneous, leakage, and protuberance drag terms were 
ignored. The empennage contribution to Co 0 was subtracted from the baseline configuration 
Cn 0 of 0.027 leaving only the wing-body parasitic drag. This allows the stabilizer areas to 
be varied and the effect on the parasitic drag to be captured. Hence, reduced stabilizer areas 
result in reduced parasitic drag. 

5.1.2 D8.2b Aerodynamic Modeling 

As a starting point, the D8.2b was modeled after the D8.1 geometry described Ref. [10]. 
The D8.2b differed from the D8.1 only in the number of engines — two engines instead of 
three — and was span limited to 118 feet, falling into the same airport operations class as 
the Boeing 737-800 [59]. An AVL model of the D8.1 concept, developed at Massachusetts 
Institute of Technology (MIT), was available from an online repository hosted at MIT. 2 An 
image of the D8.1 AVL model with visible surface loads resulting from a low-speed analysis 
can be found in Ref. [10]. 

The D8.1 model was scaled to full size (the model from MIT’s repository was 60% scale) 
so the full-scale reference areas and lengths could be used when non-dimensionalizing the 
forces and moments. The wing reference area for the D8.2b was 1110 square feet, and the lon- 
gitudinal and lateral/directional reference lengths were 11.03 feet and 118 feet respectively, 
as presented in Table 3.2. It should be noted that these reference areas were references of the 
baseline only, not actual inputs into the analysis. In the case of the multidisciplinary design 
optimization, the reference areas and lengths were recomputed with each design iteration, 
with the exception of the lateral/directional length remaining unchanged as the span was 
fixed. For the baseline AVL model used in the aerodynamic validation, the reference values 
were set according to the declared values in Fig. 3.4 for consistency. Slight discrepancies 

2 http : //web .mit . edu/drela/Publi c/web/ avl/runs/?N=A 


63 



between the declared values of Fig. 3.4 and the measured 3-view drawing values, used to 
create the OpenVSP model of Fig. 3.5, were present. 

Modeling strategies similar to the Cessna 182T AVL model were used to create the first 
D8.2b AVL model; horizontal and vertical cross sections of the fuselage were modeled as 
lifting surfaces, with the horizontal fuselage airfoils provided by the MIT AVL repository. 
The outer mold line of the fuselage, as cut by the x-y and x-z planes, was maintained by 
piecewise linear approximation using the trapezoidal sections in AVL, attempting to match 
the fuselage curvature. The D8.2b AVL model was checked to match the OpenVSP model 
as closely as possible, including all airfoils and spanwise twist. At the rear of the fuselage in 
the D8 concept, the engines were embedded to take advantage of boundary layer ingestion, 
described in Ref. [10]. Due to the complex geometry associated with the embedded engines, 
they were not included in the OpenVSP model, and were similarly ignored in the AVL model. 
Since the OpenVSP model was the outer-mold-line geometry analyzed in Cart3D, this was 
used as absolute truth in the creation of the AVL model, avoiding any unnecessary errors 
due to geometry discrepancies. 

The horizontal tail, visible in Fig. 3.5, consisted of two trapezoidal sections, one of 
constant chord between the twin vertical tails, and one of constant taper ratio external to 
the twin vertical tails, ft was desired to have simple parametric geometry for the horizontal 
tail using a single trapezoidal section. The horizontal stabilizer constant chord section was 
removed and replaced with a representative, constant taper ratio, horizontal stabilizer, and 
the single section horizontal stabilizer was used in both the OpenVSP and AVL models used 
in the aerodynamic validation. 

With the horizontal stabilizer placed near the tip of the twin verticals, care had to 
be taken to ensure the horizontal stabilizer remained connected to the verticals, especially 
during the multidisciplinary design optimization. This was accomplished by placing the 
horizontal stabilizer using the vertical tail tip leading edge, so as the vertical tail sweep 
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angle changed during the optimization the horizontal stabilizer would move as well. Consis- 
tent with transonic transports and the previous D8 study [10], the horizontal stabilizer was 
modeled as an all- moving surface for longitudinal trim in the analysis of Section 5.3.1. 

For validation of the D8 AVL model, described in Section 6.3.2, the vertical stabilizers 
were placed longitudinally and laterally using known dimensions from the OpenVSP model 
and 3- view drawing (Fig. 3.4). This laterally placed the vertical tails slightly inboard of the 
fuselage edge, and the vertical tail root chord trailing edge slightly forward of the fuselage aft 
edge. This placement was used only in the validation model for the aerodynamic comparison. 
In the AVL model used in the MDO, the trailing edge of the vertical tail root chord was 
placed at the rear edge of the fuselage to give a constant reference point, as the fuselage 
geometry remained fixed in the MDO. As the vertical root chord was increased, the vertical 
stabilizer leading edge moved forward. 

Flaps, ailerons, rudders, and an all-moving horizontal stabilizer were included as control 
surfaces in the D8 AVL model. To model the all-moving horizontal stabilizer, the hinge line 
must be placed at the leading edge, or 0% local chord, of the control surface. For a trailing 
edge device, such as a rudder, the hinge location was placed at the starting point of the control 
surface, in this case at 60% of the local chord giving a 40% local chord rudder. Ailerons were 
located on each wing with the inboard location at 67% semispan ending at 89% semispan, 
and had a control surface chord length of 20% local chord. For low-speed aerodynamics, 
flaps were modeled inboard of each aileron, starting at the edge of the fuselage and spanning 
the wing trailing edge to the start of the ailerons. This corresponds to an inboard starting 
location of 15% semispan and an outboard ending location of 67% semispan. The flap chord 
was 25% local chord. 

AVL models the control surface deflections using a theoretical effectiveness of one, mean- 
ing that any deflection results in unrealistic control effectiveness. Torenbeek shows in Ref. [60] 
that lift effectiveness for plain flaps varies nonlinearly with flap size, deflection angle, and 
gap type (open or closed). Figure 5.8 shows a carpet plot produced by Torenbeek from data 
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originating in the USAF Stability and Control Datcom [61] and Ref. [62], Olson in Ref. [63] 
discussed how the gain parameter in AVL can be used to capture the reduction in control 
surface effectiveness as a function of deflection angle. 



Figure 5.8: Lift effectiveness for plain flaps, both open and closed gap [60]. 

To accurately correlate the control surface gain with the deflection angle, especially in a 
trim analysis, a convergence loop would be required. The control surface gain parameter is 
part of the AVL input hie resulting in the aerodynamic analysis solution being insensitive to 
changes in control surface effectiveness as the surface deflections change. If the control surface 
deflection matched the gain, then the convergence could be stopped, but if the deflection 
did not match the gain, the gain was adjusted and the solution rerun until convergence was 
obtained. 

This iterative convergence was undesirable and so an alternative solution was sought. 
A fixed input gain was chosen that gave accurate control surface deflection effectiveness 
in the region where the deflection constraints would become active. For the rudders, with 
a hap chord of 40% of the local chord, a gain of 0.65 was chosen using the closed gap 
curve of Fig. 5.8. It was decided that this would be a conservative value at dehections less 
than 20 degrees, but accurate at dehections near 20 degrees, and was better than being 
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overly optimistic. The aileron gain was chosen to be 0.8. As the aileron deflections were 
envisioned to be less than the rudder deflections, a moderate gain value was selected resulting 
in pessimistic control effectiveness at angles less than ten degrees, and optimistic control 
effectiveness for deflections above 12 degrees. Finally, the flap gain was selected to be 0.75 
as only one flap angle deflection, fifteen degrees, was used in this research. 

The baseline AVL model for the D8.2b geometry is shown is Fig. 5.9 with the moment 
reference point set to 62.5 feet, as measured from the nose of the aircraft. This reference 
point was used for the baseline geometry and aerodynamic validation study only. In the 
multidisciplinary design analysis, the center of gravity, used only as a moment reference 
point in AVL, was placed a set distance, referenced to static margin, forward of the neutral 
point. This allowed the center of gravity, again the reference point as used in AVL, to shift 
with the design variables, a process described in greater detail in Section 5.3.1. In the takeoff 
flight condition, the reference point was moved 1.874 feet aft of the center of gravity and 
nine feet below the fuselage centerline, coinciding with the rear tires’ point of contact with 
the ground. 



Figure 5.9: Baseline D8.2b AVL model. 


Results discussed in Section 6.3.2 indicated that the baseline D8.2b AVL model was 
not accurately representing the aerodynamic properties as predicted in the Euler code, 
Cart3D [64,65]. The vertical cut of the fuselage remained unchanged. To better match 
the lift, drag, and pitching moment coefficients, the horizontal cut of the fuselage was re- 
moved and the wing extended to the vehicle centerline by maintaining a constant taper ratio 
and quarter-chord sweep angle. Additionally, a two degree wing incidence angle defined in 
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the OpenVSP model was removed in the AVL model to better match the Cart3D results. 
The AVL D8.2b model used for the design optimization is shown in Fig. 5.10. 



Figure 5.10: Final D8.2b AVL model used in the full system design optimization. 


5.2 Sizing and Performance Estimates using NASA’s Flight Optimization 
System (FLOPS) 

The Flight Optimization System (FLOPS) is a complete aircraft conceptual design mis- 
sion analysis and performance code that was in continuous development at NASA Langley 
Research Center starting the 1980s up until 2011 [14]. The tool is designed to take concep- 
tual level design parameters of conventional designs, such as wing area, aspect ratio, taper 
ratio, and estimate the system aerodynamics, propulsion capabilities, weights, and mission 
performance. A strength of the tool, as part of its continuous development, is its extreme 
flexibility to adjust or tune any analysis, such as adjusting wing component weight, apply 
k-factors that simulate technology factors, and optimize the mission performance with dif- 
ferent objectives and performance constraints. The high-level process flow for FLOPS is 
shown in Fig. 5.11 with the tool in analysis mode, as opposed to optimization mode, for a 
fixed design range and an externally provided engine deck. FLOPS is capable of running in 
optimization, fixed gross weight, and internally generated engine deck modes, but none of 
those features were used in this analysis. 




Figure 5.11: FLOPS process flow when run in analysis mode for a fixed design range and 
provided engine deck. 

5.2.1 Overview of FLOPS Input File 

Many of the inputs have preprogrammed default values allowing the user to successfully 
run FLOPS without explicitly including all the inputs. As more information becomes known 
about the design, such as external aerodynamics, additional inputs can be added increasing 
the size of the input hie and complexity of the model. Sectioned into namelists — FLOPS 
was written using FORTRAN 77 syntax — the input hie contains all the information required 
for running FLOPS; creating the geometry and running the analysis. The hrst namelist 
is the options namelist that controls the type of analysis being run, and hags for enabling 
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specific outputs and their filenames. Options used in this research included running the full 
analysis — contrary to just performing a weight buildup — the optimization feature turned off, 
and a detailed takeoff and landing (TOL) analysis used. 

As ordered in the input hie, the WTIN namelist is next and is where much of the 
configuration geometry information is input. Items specified in the WTIN namelist include 
the ultimate load factor, maximum landing weight in percent of the gross weight (it can 
be directly input as well), wing dihedral, vertical and horizontal tail geometry, fuselage 
information, etc. The remaining geometry information is declared in the CONFIN namelist 
where, if the optimizer feature is turned on, a range of aspect ratios, thrusts, wing areas 
and sweeps could be input. Additional information specified in the input hie includes a 
hag for the internal aerodynamic calculation or external aerodynamic override, if external 
aerodynamic data is available. Described in Section 5.3.1, both the internal aerodynamic 
analysis and external data aerodynamic override was used; the hrst FLOPS analysis used 
the internal aerodynamic analysis to compute an initial guess of the aircraft properties used 
to generate the AVL input hie, and the second analysis used the external aerodynamics 
from the cruise condition AVL aerodynamic analysis to calculate the mission performance, 
predicting total mission fuel burn. 

The payload specifications are declared in the WTIN namelist. A single-class passenger 
layout was used, affecting the furnishing weights calculated in FLOPS, with 180 passengers. 
Three flight attendants and two flight crew were added. Each passenger was assumed to 
weight 180 pounds, and have baggage weighing 45 pounds [10]. Fuselage length was specihed 
in the input, eliminating the dependency of fuselage length on passenger count. 

The MISSIN namelist is used to specify the mission structure for use in the system 
performance analysis. In the namelist, items such as taxi in/out and total takeoff time 
are input and used in FLOPS for total fuel burn calculations. Bounds on the allowable 
Mach numbers, altitude, and lift coefficient are specihed, usually by the default values, and 
the minimum fuel burn climb profile selected. Additionally, FAA limitations on indicated 


70 



airspeed are applied for altitudes less than 10,000 feet above sea level. Similar to the climb 
segment, bounds are set on the cruise segment as to maximum allowable altitude and the 
type of cruise schedule employed. A maximum allowable altitude of 41,000 feet and a fixed 
Mach number cruise, with varying altitude to optimize specific range, was used. This was the 
cruise profile that would allow for the most efficient cruise climb, reducing the total system 
fuel burn. For the descent, bounds were again placed on the allowable Mach numbers, and the 
descent profile was selected to maximize lift over drag ( L / d )• Finally, a reserve mission was 
added where 5% total fuel reserve was added to the total fuel weight, which was consistent 
with the N+3 concept D8 studies of Ref. [10]. 

5.2.2 Geometry Process 

Most parametric geometry variables are specified in the input hie, but some features, 
either not specified in the input or dependent upon multiple inputs, are calculated in the 
geometry process. Examples would be aspect ratio, if span and wing area are specified, engine 
nacelle size, multiple vertical tails (only one vertical tail geometry is input with the number 
of vertical tails as a separate input), wetted area estimates, and landing gear length. These 
calculated geometry features are passed to the aerodynamic analysis and weight estimation, 
indicated in Fig. 5.11. 

5.2.3 Aerodynamic Process 

Receiving data from both the input hie and the geometry generation process, the aero- 
dynamics used in the performance process are calculated in the aerodynamic process. The 
empirically based aerodynamic analysis, a modified version of the Empirical Drag Estima- 
tion Technique (EDET) [66], calculates the lift and drag characteristics of the input aircraft 
configuration from low speed, Mach 0.2, to the maximum input Mach number, which in 
this research was Mach 0.82. Modifications to the EDET include smoothing of the drag 
polars, increased accuracy Reynolds number calculations, and the inclusion of the Sommer 
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and Short T’ method for skin friction calculations [67]. Lift and drag polars are generated 
for a series of Mach numbers, automatically distributed by the aerodynamic analysis process, 
that are functions of the geometry and any tuning parameters that are added in the input 
hie. Examples of aerodynamic tuning capabilities include scaling lift-independent drag, lift- 
dependent drag, all subsonic drag, and all supersonic drag. Additional technology factors 
are included, such as an airfoil technology factor and the ability to include laminar how tech- 
nologies. Specihed with the laminar how technologies are the percent of local length that 
would experience laminar how, resulting in viscous drag reductions. For wings this would be 
chord, whereas for a fuselage it would be fuselage length. Laminar how was not assumed in 
this research as the D8.2b was a more current technology aircraft compared to an advanced 
technology version of the double-bubble configuration, the D8.5 [10]. 

After the lift and drag polars are generated, the information is fed to the performance 
process where any interpolations or extrapolations of the aerodynamic data are performed. 
The internal aerodynamic analysis was used in the initial FLOPS run, indicated in Fig. 5.14, 
to assist in the cruise flight condition generation and the viscous drag calculation. An option 
of the aerodynamic process in FLOPS is to override the internal aerodynamic calculations, 
and this capability was used in the second FLOPS run where AVL was used to generate the 
polars. This captured the effects of trim drag and the affect of shifting the center of gravity 
and static margin distance, but the vortex lattice code was only a potential flow analysis 
which could not capture viscous drag. In order to build up the full drag polar, the viscous 
drag calculations from the FLOPS aerodynamics analysis were used, calculated in the first 
FLOPS run, and added to the lift-dependent drag analysis of AVL. 

5.2.4 Propulsion Process 

FLOPS is capable of internally generating an engine deck using variables and control pa- 
rameters included in the input hie. Engine cycle definition decks are used for creating thrust 
and fuel how decks as functions of Mach-altitude conditions for turbojets, turboprops, mixed 
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flow turbofans, separate flow turbofans, and turbine bypass engines [67]. The FLOPS engine 
cycle analysis, described in the manual [14] , was developed and implemented by Karl Geisel- 
hart in Ref. [68]. ft was based upon the Quick Navy Engine Program (QNEP) [69], which in 
turn was based upon the Navy Engine Performance Computer Program (NEPCOMP) [70]. 

If externally generated engine data is available, FLOPS can read in the engine data 
and evaluate the mission performance. A CFM56-7BE engine deck (the engine used on the 
Boeing 737-800 aircraft 3 ) that was generated at NASA Glenn Research Center was available 
and used in this research. When the engine deck table, generated using the Numerical 
Propulsion System Simulation (NPSS) originally developed by NASA, 4 is read in by FLOPS, 
a propulsion and scaling module manipulates the data to fill in any missing data points, and 
the data is then interpolated and extrapolated, either linearly or nonlinearly, as required in 
the performance module. With the aft-mounted embedded engines of the D8 configuration, 
true freestream flow would not be ingested by the engines as typical in a pod mounted 
engine, similar to the 737-800. Modifications to the engine deck were required as the inlet 
flow to the engines included the fuselage boundary layer, called boundary layer ingestion 
(BLI). The BLI was included in the engine deck as a ram drag credit based upon the results 
from the Transport Aircraft System Optimization program (TASOPT), developed at MIT, 
and documented in Ref. [10]. 

During the system optimization, the engine static thrust was varied which required the 
scaling routine to adjust the base engine deck, giving the performance of a new engine with 
static thrust equal to the input thrust. In the final design, if the scaling was too large it 
would be best to rerun the engine cycle analysis for the adjusted thrust to obtain an engine 
model with increased accuracy. As this was only a system optimization study, no rerun 
engine decks were created for any of the optimal designs. 

3 http : //www.boeing. com 

4 http : //www. swri . org/npss/ 
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5.2.5 Weight Estimation Process 


Unlike many conceptual design tools, FLOPS does not use the traditional weight fraction 
method described by Raymer [1], Roskam [2], and Nicolai [3] for estimating gross, empty, 
and fuel weights. Geometric parameters, gross weight, aircraft type, mission range, payload, 
and maximum Mach number are all used in a set of statistical/empirical weight equations 
developed from a database of existing aircraft. The equations are dependent upon aircraft 
type, such as fighter, transport, supersonic transport, and general aviation. The transport 
aircraft empirical weight equations are based upon 17 aircraft: the C-5B, C-141B, B707-121 
B707-321, B720-022, B727-100, B737-200, B747-100, L-1011, DC-9-30, DC-10-40, B727-200, 
T-39D, L-1329, XB-70, B767-200, and B969-336C [71]. Table 5.1 shows a series of images that 
represent a majority of the aircraft used in the formulation of the transport weight equations. 5 
The weight equations were generated through a curve fit of the component weight data using 
Bayesian logic to determine the curve shape types, expanding the analysis beyond a simple 
regression [71]. Checks were included on the weight equations to ensure consistency — weight 
increases as horizontal tail size increases — and physical limits were applied — horizontal tail 
weight goes to zero as area goes to zero. Technology factors, like the aerodynamics process, 
are included to account for improved materials and manufacturing techniques compared to 
the database used to generate the weight equations. 

Data from the input hie specifies the aircraft type, in this case a transport aircraft, the 
initial gross weight guess, parametric geometry inputs, design range, payload weight and 
configuration (such as number of passengers and seat class type), and the maximum Mach 
number. Specifying the aircraft type determines which set of weight equations are used 
and the initial gross weight estimate is used to generate the system weights. For example, 
the wing weight is a function of aspect ratio, sweep, thickness, dihedral, maximum Mach 
number, and gross weight, to name a few. Using these parameters, the wing weight is 

5 All images in Table 5.1 were taken from http ://www. wikipedia.org and do not necessarily correspond 
to a specific sub-model of aircraft. For visual reference only. 
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Tabic 5.1: Reference images of aircraft used in the generation of the FLOPS transport weight 
equations. 



B727 





DC-10 


XB-70 




B767 



B720 




L-1329 


estimated, including the primary and secondary structure. Engines mounted on the wings 
require additional wing weight but also provide inertia relief, both taken into account in the 
weight equations. The number of passengers, their class configuration, and weight is used 
in the fuselage weight equation, consisting of structure and passenger amenities. A small 
sensitivity to range is included in the weight equations as long range missions have increased 
passenger amenities compared to short range missions, even for an identical aircraft type. 

Weight groups, as output in FLOPS, were a combination of system aircraft components 
where the groups were combined using typical design level categories. Table 5.2 summarizes 
the output weight statement from FLOPS [68]. 
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Table 5.2: Aircraft weight statement summary as output from FLOPS [68]. 

Group 

Components 

Structure 

Propulsion 

System & Equipment 

Operational 

Payload 

Wing, HT, VT, Fuselage, Landing Gear, Nacelle 
Engines, Thrust Reversers, Misc. Systems, Fuel System 
Surface Controls, APU, Instruments, Hydraulics, Electri- 
cal, Avionics, Furnishings, AC, Anti-icing 
Crew & Baggage, Unusable Fuel, Engine Oil, Passenger 
Service, Cargo Containers 
Passengers, Baggage, Payload 


5.2.6 Mission Analysis Process 


Composed of multiple mission segments specified by the user, the configuration mission 
performance is evaluated through a weight-based integration where the total mission fuel 
burn, range, and detailed flight profiles are calculated. The mission profiles are generated 


by a sequential input of segments, up to a maximum of 40 segments, with the segment input 
options indicated in Table 5.3. The analysis must always begin with a start segment, finish 
with an end segment, and contain at least one cruise segment. 

Table 5.3: Available mission segments in FLOPS along with their primary inputs [68]. 

Segment 

Primary Input 

Start 

Starting Mach number and altitude 

Climb 

Climb schedule number 

Cruise 

Cruise schedule number and total distance to end 

Refuel 

Fuel added and time required 

Release 

Weight released 

Accel 

Engine power setting and ending Mach number 

Turn 

Turn arc and engine power setting or turn acceleration 

Hold 

Cruise schedule number and time 

Descent 

Descent schedule 

End 

Ending Mach number and altitude 


Figure 5.12 shows two example mission profiles that would typically be analyzed in 
FLOPS, and two key features worth noting. They are the marked “free segment” and the 
“instantaneous descent” . Each mission profile must include a cruise segment that does not 
have a specified segment range. FLOPS does not “fly” the mission from beginning to end. 
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Instead, the analysis starts with aircraft ramp weight, subtracts taxi out and takeoff fuel 
allowances, and then begins stepping through the mission segments until the free cruise 
segment is reached. The analysis then skips to the end of the profile and starts with the zero 
fuel weight, adds the reserve fuel, and then steps through the descent, cruise, hold, and other 
segments in reverse until the free cruise segment is reached from the opposite direction. The 
difference in fuel between both sides of the free segment then determines the “free” cruise 
segment distance. 



Figure 5.12: Example mission profiles as analyzed in FLOPS [68]. 


Several options are available to optimize the climb profile, depending on the desired 
mission performance. Allowable profiles include minimum fuel- and time-to-distance profiles, 
or a weighted combination of both. Instead of distance, minimum time- and fuel-to-climb can 
also be selected, again with the option of having a weighted combination of the two. Table 5.4 
summarizes the different climb options. Applications of the different optimal climb profiles 
are suggested by the FLOPS manual, particularly to which type of vehicle configuration the 
climb profile should be applied. Climb profiles optimized on minimum fuel-to-distance are 
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applicable to subsonic transports but should not be used for supersonic transports; minimum 
fuel-to-climb should only be used for interceptors; minimum time-to-climb is suggested only 
for fighters; both subsonic and supersonic transports can be optimized on minimum fuel- 
to-climb [67]. Any combination of climb profiles can be used in a mission profile, up to a 
maximum of ten climb segments with FAA climb restrictions below 10,000 feet applied. 

Table 5.4: Climb schedule options for use in FLOPS climb segment analysis [67]. 


Flag 

Description 

1 

Minimum fuel-to-distance 

0 

Minimum time-to-distance 

0-1 

Combination of minimum fuel- and distance-to-climb 

-0.001 

Minimum time-to-climb 

-1 

Minimum fuel-to-climb 

- 0 . 001 - -1 

Combination of minimum time- and fuel-to-climb 


Cruise schedules can be optimized with more options than both the climb and descent 
schedules, with the options summarized in Table 5.5. Unlike the climb schedule, cruise 
schedule optimizations are discrete and a weighted combination of two options is not allowed. 
The cruise schedule options include different possible combinations of fixed Mach number, 
fixed altitude, optimal Mach number, optimal altitude, and maximum Mach number, which 
are optimized for either segment endurance or specific range, velocity divided by fuel flow. 
A fixed lift coefficient cruise schedule can also be selected. If a hold segment is specified, 
the cruise schedule options are applied to the hold segment. A maximum of fifteen cruise 
segments can be used in a mission profile, again any combination of cruise schedules can be 
chosen. Reserve cruise segments are included in the cruise schedule definitions [67, 68] . 

Descent segments in the mission profile have three schedule options, summarized in 
Table 5.6. In the case of military aircraft, no time, distance, or fuel credit is given for 
descent segments, so an option for no-credit descent is available with a zero flag. For all 
other aircraft types, a descent at maximum lift-to-drag ratio or constant lift coefficient can 
be selected, and FAA restrictions on calibrated airspeed below 10,000 feet applied, same 
as the climb schedule. Limits in the maximum dynamic pressure can also be chosen for 
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Table 5.5: Cruise schedule options for use in FLOPS cruise segment analysis [67]. 
Flag Description 

0 Optimum altitude and Mach number for specific range 

1 Fixed Mach number and optimal altitude for specific range 

2 Fixed Mach number at input maximum altitude or cruise ceiling 

3 Fixed altitude and optimum Mach number for specific range 

4 Fixed altitude and optimum Mach number for endurance 

5 Fixed altitude at a constant lift coefficient 

6 Fixed Mach number and optimum altitude for endurance 

7 Optimum Mach number and altitude for endurance 

8 Maximum Mach number at input fixed altitude 

9 Maximum Mach number at optimum altitude 

10 Fixed Mach number at constant lift coefficient 


the descent schedule. A special property of the descent segment only, any descent segment 
followed by a climb segment will be considered, for all but one of those descent segments, an 
instantaneous (zero fuel, time, and distance) change in altitude [68]. 

Table 5.6: Descent schedule options for use in FLOPS descent segment analysis [67]. 
Flag Description 

0 No descent time, distance, or fuel credits 

1 Descent at optimum lift-to-drag ratio 

2 Descent at constant lift coefficient 


With the primary mission completed, additional fuel must be added to account for 
operation variations due to weather, missed approach, etc. How the reserve fuel is calculated 
depends upon the type of reserve mission used, with various options available in FLOPS. 
Two variables, RESRFU and RESTRP, are used in the calculation of constant fuel reserves. 
If RESRFU > 1, a fixed reserve fuel weight is added to the mission fuel weight where 
the value is input in pounds, and if RESRFU < 1, a fraction of total usable fuel weight 
was added. RESTRP is used to add reserve fuel as a fraction of total trip fuel weight. An 
alternate airport can be used and the distance to the alternate input with the ALTRAN 
variable. Combinations of RESRFU, RESTRP, and ALTRAN can be used in the definition 
of the complete mission reserve [67]. 
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Reserve fuel can be calculated one of three ways: 1) calculate reserve fuel required for 
trip to alternate airport plus RESRFU and/or RESTRP, 2) reserve fuel calculated using 
RESRFU and RESTRP only, and 3) the reserve fuel is the remainder after the primary 
mission has been flown [67]. Additional fuel is added for an input missed approach time. 

This research used a typical transport aircraft mission profile with segments consisting 
of climb, cruise, and descent. The mission segments used in the FLOPS input hie were 
START, CLIMB, CRUISE, DESCENT, and END, with a fuel reserve added using the RE- 
STRP variable set to 5%. The single cruise segment was defined as a free segment. This 
profile matched the mission profile of the N+3 advanced concept study, defining the D8 
configuration, from Ref. [10]. As minimization of total fuel burn was the objective of this 
research, minimum fuel-to-climb was the chosen schedule for the climb segment, restricted 
to the FAA operational constraint limiting calibrated airspeed to 250 knots for altitudes less 
than 10,000 feet. Cruise Mach number was specified in Ref. [10] and was used again here, 
set to Mach 0.72. With the cruise Mach number fixed, and an objective of minimizing fuel 
burn, the fixed Mach number with optimum altitude for specific range was the selected cruise 
schedule, allowing for a cruise climb. Cruise altitude was limited to a maximum of 41,000 
feet. Maximum lift-to-drag ratio was selected for the descent schedule, with a two minute 
allowance for a missed approach. 

5.2.7 Balanced Mission Analysis 

Recalling Fig. 5.14, based upon an initial guess for gross weight the aircraft empty weight 
is calculated, and the difference between the operational weight fully loaded and gross weight 
is assumed to be fuel. The mission is then flown with the free segment flown with only the 
fuel remaining after all other mission segments, determining total range. In the case of a 
design range mission, if the calculated range is equal to the design range, the analysis moves 
on to the takeoff and landing analysis, a completely separate analysis process in FLOPS. Any 
difference between the design range and the calculated range from the performance process, 
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the gross weight estimate will be updated and the analysis run again until convergence is 
achieved on the design mission range. All mission fuel, less defined reserves and taxi fuel, 
will be burned, while achieving the design range, leaving zero mission fuel at the end of the 
mission. 

5.2.8 Takeoff and Landing Analysis 

The detailed takeoff and landing analysis procedure was used in the multidisciplinary 
design optimization. Developed as a physics-based, first principles analysis, all FAR Part 
25 (civilian transports) or MIL-STD-1793 (military aircraft) airworthiness requirements are 
applied. During takeoff, the time- integrated analysis captures the variation in thrust with 
velocity and the changes in lift and drag coefficients as the aircraft rotates, and lifts off. 
During the descent the aircraft will flare, which changes the lift and drag coefficients, touch 
down, deploy spoilers if included, and apply the brakes. In the takeoff analysis, numerous 
simulations are run in order to determine the balanced held length which is the greatest 
distance of the following: 

• one engine inoperative (OEI) held length to clear 35 foot obstacle 

• aborted takeoff at decision speed with one engine inoperative 

• aborted takeoff at decision speed with all engines operating (AEO) 

• 115% of distance to clear 35 foot obstacle with all engines operating 

Excess thrust at the second segment climb gradient, dependent upon the number of engines 
in the event of an engine failure, must be greater than zero in order to have a successful 
takeoff. 

The maximum takeoff distance allowed in the multidisciplinary design optimization 
was set to 8,000 feet. This was chosen to remain consistent with the previous research on 
the D8 geometry prior to the application of a 5,000 feet balanced held length requirement 
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discussed in Ref. [10]. For the landing field length requirement, the Boeing Commercial 
Aircraft Company document for airport planning entitled, “737: Airplane Characteristics 
for Airport Planning” was used [72], Standard day, sea level conditions were used at the 
maximum allowable landing weight, 146,300 pounds. Figure 5.13 shows at the maximum 
landing weight, the held length for a dry runway at sea level, standard day is 5,800 feet. This 
was chosen as the landing field requirement for the multidisciplinary design optimization. 
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Figure 5.13: Boeing 737-800 landing held length plot for both wet and dry fields, various 
altitudes, at standard day temperature [72], 
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5.3 D8.2b System Analysis and Multidisciplinary Design Optimization 

Model Center® , a tool used within the Aeronautics Systems Analysis Branch at NASA 
Langley Research Center, was used to integrate the numerous disciplines that have been dis- 
cussed to this point. Developed by Phoenix Integration, 6 ModelCenter is a multidisciplinary 
integration, design, and optimization software that allows the user to integrate their various 
analysis tools methods while facilitating all data transfer. The design variables are speci- 
fied and can be brought into any of the integrated optimization schemes, including gradient 
methods, response surface optimizers, genetic algorithms, and particle swarm optimizers. 
During an optimization, multiple instances are able to be run allowing for the parallelization 
of the analysis (if the optimization scheme allows), drastically decreasing optimization run 
time. 

5.3.1 Integrated Analysis in ModelCenter® 

FLOPS, the AVL input hie generation, AVL, the flight condition calculations, and the 
MATLAB® dynamic analysis have all been seamlessly integrated into ModelCenter® to 
perform the multidisciplinary design analysis and optimization. Data was sent from each 
analysis component to the next through the linking capability in ModelCenter, where the 
output of a component was linked to the inputs of one or more components. This allowed 
components to be run in parallel, accelerating the analysis process. Additionally, using the 
Analysis Server software in collaboration with ModelCenter, multiple instances of the model 
were run allowing for multiple designs to be evaluated simultaneously, reducing the total 
time for the complex design optimization. 

Figure 5.14 depicts the multidisciplinary design analysis, as implemented in Model- 
Center, in a process diagram. The first step was the initialization of the design variables, 
controlled either by the user or by the optimization routine. The design variables allowed to 
vary were wing area, wing quarter-chord sweep, dihedral, wing longitudinal apex location, 

6 http : / / www . phoenix- int . com/ 
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static thrust, horizontal tail area, vertical tail area, and vertical tail sweep. An additional 
design variable that was allowed to vary was the minimum static margin, depending on the 
design case. For the cases where minimum static margin was not a variable, it was hxed at 
10% of the mean aerodynamic chord. 

No weight buildup was used for the calculation of center of gravity. To allow the 
optimizer the flexibility to place the center of gravity anywhere on the configuration, the CG 
was placed based upon the aircraft neutral point and specified static margin. The aircraft 
neutral point was a function of the design variables allowing the optimizer to shift the neutral 
point, and by association, the center of gravity placement. This was chosen to open the 
design space, freeing the optimizer to place the center of gravity as desired, and eliminating 
the guesswork in placing each component center of gravity, including subsystems. This 
allows for potentially unique designs that would not be eliminated dne to center of gravity 
constraints. As advanced air transport technologies evolve, high mass density subsystems, 
such as batteries, can be placed to allow for a larger range of center of gravity placements 
for a particular configuration. A summary of the design variables, and their allowable range, 
is summarized in Table 5.7. The ranges were selected as to not become active and corner 
the optimizer in the design space. 


Table 5.7: Summary of design variables allowed to be varied by the optimizer with their 
allowable ranges. Static margin was only allowed to vary in select design cases. 


Design Variable 

Description 

Nominal 

Minimum 

Maximum 

T 

Thrust (lb.) 

25,000 

10,000 

30,000 

S 

Wing Area (sq. ft) 

1100 

800 

2500 

A 

Wing Sweep (deg) 

10 

0 

30 

r 

Wing Dihedral (deg) 

5 

0 

10 

Sht 

HT Area (sq. ft) 

277 

100 

500 

SvT 

VT Area (sq. ft) 

74 

50 

300 

A V t 

VT Sweep (deg) 

40 

0 

65 

X w 

apex 

Wing Apex (ft) 

55 

20 

75 

SM 

Static Margin 

10% 

-5% 

30% 
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Figure 5.14: Multidisciplinary design analysis process chart with data flow direction indicated 
by the arrows. 
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For the horizontal stabilizer, only the planform area was allowed to be a design variable. 
Located at the tip of the twin vertical stabilizers, allowing the horizontal tail sweep to be 
a variable would disconnect the horizontal stabilizer from the tips of the verticals, resulting 
in a stabilizer floating in space. To ensure this did not happen, the design variables for 
the horizontal stabilizer were restricted, with any positional movement implemented by the 
vertical stabilizer geometry. 

Wing aspect ratio was not allowed to be a design variable due to the fixed span re- 
quirement of the D8.2b configuration. The goal of the D8.2b configuration was to quantify 
some of the configuration and technology impacts while maintaining similar operational ca- 
pabilities as the Boeing 737-800 aircraft. Increases in span would provide large induced drag 
benefits while moving the vehicle to a larger airport operations class. Also, increasing the 
airport operational class would cause the configuration to be non-comparable to the 737-800 
in terms of operational destinations. Additionally, as aspect ratio, and therefore span, is 
increased, stress-based constraints become inactive in the design and are replaced by flutter 
constraints, an analysis well beyond the scope of this research. Therefore, the span was fixed 
at 118 feet, same as the 737-800, to eliminate these issues. 

Similarly to the wing, the horizontal tail aspect ratio was also kept constant. This was 
done to reduce the number of design variables, limit the span due to flutter concerns, and 
maintain a constant lift curve slope as described in Section 4.5.3. 

ModelCenter, through the use of an custom component called a Script Wrapper, created 
the FLOPS input hie and executed FLOPS. This first FLOPS execution was used to extract 
dependent geometry parameters that were required by AVL. To ensure consistency in the 
meta-geometries between all components, all the design parameters were linked to eliminate 
the possibility of failing to update an input in a later analysis. In addition to generating 
geometry information required for the AVL input hie, FLOPS was used to estimate the sys- 
tem weights, takeoff rotation velocity, and start of cruise flight conditions. This information 



was used in the creation of the different flight conditions for the static and dynamic analysis, 
described in Section 4.4.2. 

Any inputs not required by FLOPS, but necessary in the generation of an AVL input 
hie, were directly sent to the AVL Geometry Generation component. Again, a custom 
Script Wrapper created for this research, the AVL geometry generation component created 
the input hie used in all AVL analyses. Required parameters for the generation of the AVL 
input hie included the control surface sizes and locations, and the longitudinal wing root, 
leading edge apex location, X Wapex . 

The initial AVL analysis was used for two reasons: 1) to capture the lift coefficient and 
zero geometric angle of attack with dehected haps in ground effect, and 2) to calculate the 
aircraft neutral point for placement of the center of gravity based upon the static margin 
design variable. Instead of using a hxed, approximated lift coefficient for the takeoff analysis, 
the initial run of AVL was able to predict the lift coefficient in ground effect, taking into 
account the change in lift due to the hap deflection. This made the lift coefficient used in 
the takeoff flight condition unique to each configuration in the design optimization. Also, 
unique to this research was the use of static margin as a design variable, where the center of 
gravity was placed based upon the neutral point and specihed static margin instead of the 
typical calculation of center of gravity. To place the center of gravity, an initial run of the 
AVL model had to be used to determine the configuration neutral point. 

The standard atmosphere component was used in the generation of the hve flight con- 
ditions described in Section 4.4.2. Temperature, pressure, and density were calculated, as 
functions of altitude, using the standard atmosphere equations presented by Anderson [73]. 
Viscosity was calculated using Sutherland’s formula. 7 With known altitude at each flight 
condition, determined from the initial FLOPS run, and calculated temperature, the speed 
of sound was calculated. 

' http : / / www . gre . nasa . gov/WWW/k- 12/airplane/ viscosity . html 
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The five flight conditions, discussed in Section 4.4.2, were calculated in the flight condi- 
tions component. Data calculated in this component were passed to the full AVL analysis for 
evaluation of the static and dynamic vehicle performance, using steady state trim deflections 
and stability derivatives. Placed off the specified static margin, either fixed or as an opti- 
mization design variable, the rearmost center of gravity position was calculated in the CG 
placement component and fed into the AVL analysis and FLOPS aerodynamic override. The 
main landing gear position, placed a static distance of 1.874 feet aft of the rearmost center 
of gravity position and nine feet below the aircraft centerline (Section 5.1.2), was calculated 
in the takeoff flight condition and placed in the CG/Landing Gear Placement component. 
This set the reference point location for the takeoff trim analysis in AVL. 

In the AVL analysis component, five flight conditions were analyzed with the moments 
trimmed to zero, resulting in a steady-state condition. The exception was the takeoff analysis 
pitching moment, which was trimmed to the required aerodynamic pitching moment for 
lifting the nose wheel off the ground, calculated in the takeoff flight condition. In the static 
flight conditions — takeoff, one engine inoperative, maneuver — the trimmed control surface 
deflections were output for evaluation against the static constraints, which were used in the 
evaluation of the objective function. In the dynamic flight conditions of stall and cruise, 
the steady-state angle of attack control surface deflections were output into the dynamic 
analysis. Additionally, the stability derivatives were calculated in AVL and passed to the 
dynamic analysis for use in the linearized perturbation equations of motion. 

For both the stall and cruise flight conditions, the calculated stability derivatives, steady- 
state angle of attack, and trim deflections were input into the system dynamic equations of 
Section 4.1 in the Dynamic Analysis component. The optimal control gains were calculated 
and the system perturbation, atmospheric turbulence, and discrete gust responses were all 
evaluated in the dynamic analysis, with the detailed methodology described in Chapter 4. 
The stall and cruise flight conditions were analyzed independently allowing each flight condi- 
tion to be evaluated in isolation against the dynamic constraint requirements. Violations of 



the constraints from the dynamic analysis were included in the objective function as penalties 
on the objective function, or otherwise fuel burn. 

FLOPS contained a procedure for estimating the aerodynamic performance of the input 
parametric geometry, discussed in detail in Section 5.2. Given a fixed geometry, adjusting 
the static margin and center of gravity positions has no effect on the aerodynamic analysis in 
FLOPS as the analysis was insensitive to variations in drag due to trim or angle of attack. To 
capture the sensitivity of induced and trim drag due to changes in geometry, static margin, 
and center of gravity position, AVL was again used to create a series of drag polars, as a 
function of Mach number, to be included in the FLOPS input hie as an aerodynamic override. 
A sweep of lift coefficients ranging from 0.0 to 0.6 in increments of 0.2 were used with three 
Mach numbers: a low-speed Mach number of 0.4, a mid-speed Mach number of 0.55, and 
the cruise Mach number of 0.72. Overriding the internal aerodynamic calculation with the 
aerodynamic calculation of AVL captured the configuration drag coefficient sensitivity to 
static margin and center of gravity location. 

With the aerodynamics calculated in AVL, FLOPS was run a second time with the 
same input hie as the hrst FLOPS run at the beginning of the analysis process, with the 
exception of the external aerodynamics substituted for the internal analysis. The second 
FLOPS analysis was used to capture the changes in fuel burn due to changes in trim drag, 
with the fuel burn output used as the minimized objective in the system optimization. 

The overarching objective was to minimize system fuel burn, which was output from 
the second FLOPS analysis. Any violations of the static or dynamic constraints, or any 
performance constraints violations, resulted in a fuel burn penalty. Moderate violations of 
the constraints were possible in the case where the objective function, with the constraint 
violation penalty, would be less than the objective function with no penalties. 



5.3.2 Optimization Methodology 


Many different optimization routines and methodologies exist for determining an “opti- 
mal” design, such as gradient-based optimizers, surrogate model optimizers (such as Boeing’s 
proprietary optimizer, Design Explorer, 8 ) particle swarm optimizers, and genetic algorithm 
optimizers. The challenge lies in guaranteeing that the found optimum is the global opti- 
mum, not a local optimum in a design space encompassing multiple minima. Many times 
this guarantee is not possible. Thus, the goal of any optimization is to find a feasible design 
that minimizes the objective function while any improvements from the found design would 
be minimal. 

Each optimization routine has benefits and detriments resulting in a trade-off between 
convergence speed and robustness of the routine in finding the optimal solution. Gradient- 
based optimizers have the benefit of converging to the optimal solution in a smaller number 
of iterations compared to other routines. The gradient-based optimizer is ideal in smooth 
and continuous design spaces, calculating a mathematically-provable global optimal design 
in a concave design space. However, a realistic design space is rarely concave and any 
design space encompassing functions greater than second order, with changes in concavity, 
will introduce local minima. Some optimization routines are driven to these local minima, 
dependent upon the sensitivity to initial conditions and the value of the local gradient. 
Additionally, gradient-based routines are not applicable to a discontinuous design space. 

Generally, the next iteration in a gradient-based optimization routine depends upon 
the solution of the previous iteration, which does not lend itself well to parallel processing. 
If the solution execution time is small, this is not a problem, but if the solution time is 
large, then this can be quite limiting. The lack of parallel computing capability is made 
up by the speed of convergence, but only if the derivatives of the design space are available 
with respect to the design variables. If not, a method of computing the derivatives, such as 
finite difference, must be used which can result in an excess number of analysis executions, 
8 http : //www.boeing. com/ news /frontiers /archive/2005/ april/ i_tt . html 
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slowing the optimization process. Typically, the gradient-based optimizer is used when the 
execution time of the analysis is long and the design space is well-behaved, meaning it is 
smooth, continuous and local minima only differ slightly from the global minimum. 

A problem with gradient-based optimizations in ModelCenter® is that the ModelCenter 
routines are incapable of handling analysis failures. Instances of the analysis may fail, due 
to the optimizer exploring extreme regions of the design space. The failure results in data 
output hies that cannot be successfully parsed, extracting the required data to continue the 
optimization process. As an example, a case could be run in AVL that would be unable 
to trim, resulting from a trim convergence failure, and results of the analysis would not be 
parsed correctly. 

A robust optimizer that can handle analysis failures, such as a surrogate or genetic 
algorithm optimizer, is required. It is desired to maintain fast convergence rates, while 
reducing the number of iterations to convergence. Boeing’s Design Explorer, a surrogate 
model optimizer, is ideal for such an analysis. Design Explorer is designed to efficiently solve 
engineering problems where the analysis code run times arc long, the design space is noisy and 
non-smooth, and failures can occur in the analysis. Additionally, the algorithm is designed 
to be less likely to stop at a local minima resulting in larger success rates in Ending the global 
minimum [74], Figure 5.15 shows the process how of the Design Explorer optimizers. The 
process is started by carefully choosing values of the design variables, exploring the design 
space using an orthogonal array [74], then the analysis code is run to begin to develop a 
response model capturing the relevant responses. 

From the results of the initial runs, surrogate models are developed to approximate the 
analysis with arbitrary values of the design variables, within the limits of the lower and 
upper bounds. These surrogate models are designed to be evaluated quickly, in comparison 
to the analysis code, and are continuous in order to predict the value of the objective as 
if the analysis code was run. The surrogate models are created using interpolating Kriging 
models, where they are able to be recalibrated as new data becomes available [74-76]. With 
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Figure 5.15: Process flow for Design Explorer optimizer as shown in the ModelCenter help 
menu [74], 


the generated surrogate models, a gradient-based optimizer is run using the surrogate models 
numerous times from randomly generated starting points. 

The results of the surrogate model optimization are one or more unique local minima. 
These solutions are compared to the solution of the analysis code run with the same design 
variable values. For simple problems, this comparison can be quite good after very few iter- 
ations. For complex design problems, the error between the analysis code and the surrogate 
model optimization can be large. The actual values obtained from running the analysis are 
used to tune, or calibrate, the surrogate model and the gradient-based optimizer is rerun. 
This process is repeated until either the maximums of the optimizer settings are reached, 
as input by the user (such as maximum function evaluations), or no improved designs have 
been found [74], 

When no improved designs have been found from the iterative process of refining the 
surrogate model and running the gradient optimizer, a local pattern search is performed to 
explore the design space in the vicinity of the found optimum. If no improved designs are 
found, then the optimization process is terminated and the design is guaranteed to be at 
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least a local minim um. If a better solution is found, the local pattern search is terminated 
and the surrogate model refinement process is restarted [74], 

Design Explorer was the preferred optimizer for this research as the execution time of 
the full multidisciplinary analysis was long. Additionally, Design Explorer is capable of being 
run with parallel processors, thus taking advantage of ModelCenter’s capability of running 
multiple instances of the model, drastically reducing overall computation time. Unfortu- 
nately, the design space was extremely noisy, discontinuous, and non-smooth, preventing the 
Design Explorer algorithm from generating surrogate models that accurately mapped the 
design space. As a result, the Design Explorer optimizer could not be used in the system 
optimization. It is possible that Design Explorer could be expertly set up to successfully 
perform the optimization with different algorithm settings, but the author was unsuccess- 
ful. As a result, computation time was sacrificed for algorithm robustness with the use of a 
genetic algorithm optimizer. 

Optimization using the Darwin Algorithm 

Darwin, a genetic algorithm search routine developed by Phoenix Integration, Inc., 9 
was the optimizer used in the multidisciplinary design optimization of this research. The 
algorithm is capable of using both discrete and continuous design variables with any specified 
number of constraints. Single objective optimization problems reach convergence when the 
fitness function has not improved over a specified number of generations, whereas a multi- 
objective optimization problem uses a Pareto front in determining convergence [74,77]. The 
minimization of fuel burn was used as a single objective optimization problem. To determine 
convergence, a specified number of sequential generations without objective improvement was 
selected with a limit on the number of total generations to be analyzed. 

The Darwin algorithm control options dialog box, as integrated in ModelCenter®, is 
shown in Fig. 5.16. Many of the options were left at the system defaults, such as population 

9 http : / / phoenix- int . com 
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size. The population size, a function of the number of design variables, was set to 63 
and was constant throughout the optimization. Decreasing the population size reduces the 
optimization run time at the risk of being caught in local minima, a result of a lack of design 
diversity in the population [74], Increasing the population size will increase the number of 
iterations, but will be more robust in avoiding local minima. The default value for population 
size, based on the number of design variables, was used in all optimization cases. 



Figure 5.16: Darwin algorithm options as captured in a screen shot from ModelCenter. 
Options shown are the values used in the full static trim and dynamic optimization cases. 


Darwin has the option of two different selection schemes: elitist and multiple elitist. 
Multiple elitist was chosen as it is more effective in problems where many local minima 
surround the global minima. Since the specific topography of the design space was unknown, 
this was deemed more robust. The multiple elitist selection scheme works by combining the 
parent and child populations into one list and ranking them based upon their fitness, the 
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sum of the objective function and any penalties [74], The next generation is created by 
taking the top N p number of designs, where N p is the specified number of preserved designs, 
to the new generation and filling the remaining designs with the top ranked child designs 
that have not already been selected. The number of preserved designs, N p , should be kept 
small as that introduces the greatest number of new designs. Large values of N p reduce the 
number of new designs and can result in the algorithm search becoming localized, trapping 
the search in a local minimum [74], 

The initial population was generated randomly through the use of a seed value that 
was set to zero, which allows ModelCenter to randomly generate the seed. This results in 
successive runs not necessarily giving the same results, or if the result was the same, taking 
a different path to get there. However, if it is desired to recreate an optimization run, 
the randomly generated seed value is saved as part of the output during the optimization 
run. Inputting the same seed value, as long as all other settings remain unchanged, will 
reproduce the same optimization run. Darwin’s memory function was also used with the 
goal of improving the efficiency of the design optimization [74], A discrete design, along 
with its response, is stored in a binary tree, eliminating the need to reanalyze duplicate 
designs that are discovered in the optimization process [74], 

The convergence criteria, as mentioned previously, differ between single objective and 
multi-objective optimizations. A single objective can be minimized, creating a best design, 
whereas a multi-objective optimization has no best design, but rather a trade-off between 
the different objectives, typically shown as a Pareto front. The Darwin algorithm can be 
stopped in two ways: stop after a fixed number of generations, or stop after the solution 
has converged without any improvement for a specified number of generations. The second 
option was used as the convergence criteria, where the number of sequential generations 
without improvement ranged from 10-20, depending on the run case; the larger the number 
of design variables and constraints, the larger the number of successive generations without 
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improvement required for convergence. The maximum number of generations was reduced 
to 200 from the default of 1,000 to bound the runtime. 

Along with the multiple elitist selection scheme, crossover and mutation are used in 
the creation of the next generation. A uniform crossover procedure was applied with a high 
probability, typically 0.8 < P c < 1.0, as suggested by the Darwin algorithm manual and 
ModelCenter default value, to traverse the design space [74,77]. The child designs subjected 
to crossover are forced to be different than all other child and parent designs, and the process 
is repeated to fill the population of the next generation. Mutation of the design’s genetic 
string introduces random alterations into the population while preventing premature loss 
of important genetic information [74,77]. It also brings in design features that may have 
never been represented by the initial population, diversifying the overall design population. 
During mutation, a single value in the genetic string representing a particular design is 
changed, at random, to any other permissible value. This mutation process, applied with a 
low probability of 0.01 < P m < 0.3, occurs after the crossover operation and completes the 
creation of a new generation in the optimization process [74,77], summarized by the process 
chart of Fig. 5.17. 

Constraint tolerance is the last set of options for the Darwin algorithm displayed in the 
options toolbox window of Fig. 5.16. The designs in each generation, and each successive 
generation, are ranked according to their fitness. Calculated by Eq. 5.3 [74], the fitness, /, 
is a function of the objective function, o, and any penalties p due to violated constraints. 

f = o + p (5.3) 

The “maximum constraint violation” and “percent penalty” options in the dialog box give 
the user control over how the penalty function is applied. Described in Eq. 5.4 [74], the 
penalty function only provides a slight penalty to violated constraints where the violation is 
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Figure 5.17: Darwin genetic algorithm process flow chart, recreated from Ref. [77]. 
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less than the user specified maximum, in percent, 


p = p* 



(5.4) 


where p* is the percent penalty, e t is the total constraint violation, and e m is the maximum 
allowed constraint violation. In the case of Fig. 5.16, any constraint violation less than 5% 
will be only lightly penalized, but as the violation increases beyond the maximum allowed, the 
penalty increases drastically. If increased margin on the constraints is desired, the percent 
penalty weighting (50% in Fig. 5.16) should be reduced and/or the maximum constraint 
violation increased. To the contrary, the opposite changes would be made to tighten the 
tolerances on the optimization constraints. The default constraint tolerance values, 50% 
penalty, and 5% maximum constraint violation, were used in all optimization cases of this 


research. 



Chapter 6 

Verification and/or Validation of Methodology and Analysis Tools 

6.1 Verification of Equations of Motion 

The fully coupled nonlinear equations of motion were derived with the aide of the 
MATLAB® Symbolic Toolbox with both the aerodynamic and moment terms approximated 
using a Taylor series expansion, dropping all terms over first order. The derived perturbation 
equations were verified against Roskam, Napolitano, and Schmidt in Refs. [15,34,78]. 

An eigenanalysis was performed on the state matrix of the derived equations of motion 
to compare roots of the characteristic equation given in Napolitano [15] and eigenvalues cal- 
culated internally in AVL. The stability derivatives given in Refs. [15,34] were input into the 
full-coupled perturbation equations and the eigenvalues of the state matrix were computed. 
The eigenvalues are compared in Fig. 6.1. The AVL calculated stability derivatives were 
input into the perturbation equations and the eigenvalues calculated for both methods. The 
results are plotted in Fig. 6.2. In both cases, the eigenvalues match with the input sta- 
bility derivatives coming from separate sources indicating the derived equations of motion 
accurately represent the vehicle dynamics. 
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Figure 6.1: Derived equations of motion modes compared to modes presented by Napoli- 
tano [15]. 



Figure 6.2: AVL computed modes compared modes of derived equations of motion. 
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6.2 Verification of the Atmospheric Disturbances 


The von Karman continuous turbulence spectrum was calculated for varying length 
scales to ensure the accurate calculation of the spectrnms. Length scales of 500, 1,000, 2,000, 
and 5,000 feet were used in the vertical spectrum equation for varying spatial frequencies, 
ff, between 10 -4 to 10” 1 . Figure 6.3 shows the vertical spectra calculated for the different 
length scales and these results agree with Fig. 9.56 in Ref. [30]. 



Figure 6.3: The von Karman continuous turbulence spectrum. 

LIsing the vertical spectra of Fig. 6.3, the results of the example in Section 13.4 of 
Ref. [29] were accurately recreated for the set of equations of motion, and appropriate sta- 
bility derivatives, given in the example. In the example, the point approximation, applying 
the gust to the center of gravity of the aircraft, works very well in capturing the response of 
the motion states, a , 6 , and u. 
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6.3 Validation of the Aerodynamic Analysis 


Validation of the aerodynamic model was necessary to ensure the stability derivatives 
used in the dynamic analysis were a reasonable reflection of reality. In conceptual design 
where the design space is being rapidly explored, it is unrealistic to use high-order tools 
(CFD) and/or experimental testing to formulate the aerodynamic model. With each change 
in the geometry new aerodynamic models have to be created; the lack of fluidity in high- 
order tools and experimental testing eliminates this as a possibility with current capabilities. 
Lower-order methods must be used that can rapidly estimate the aerodynamic properties 
of a geometry while matching high-order tools or experimental results reasonably well. An 
accurate representation of the geometry is desired while keeping in mind the limitations of 
the low-order methods. 

Experimental stability derivatives of the Cessna 182 in the cruise condition were avail- 
able through Ref. [34] and were used to validate the Cessna 182T AVL model. The D8.2b 
geometry, due to the conceptual nature of the geometry, has no experimental data for com- 
parison and so the AVL model was compared to results of Cart3D, a high-fidelity inviscid 
analysis CFD tool using Cartesian mesh methods developed jointly by NASA Ames and 
the Courant Institute at NYU. 1 The Cart3D analysis was run at multiple Mach numbers to 
capture the aerodynamic performance throughout the climb, cruise, decent mission segments. 

6.3.1 Cessna 182T Aerodynamic Model 

Stability derivative data at multiple flight conditions for the Cessna 182T were used 
for validation of the AVL model described in Section 5.1. As previously described, the AVL 
geometry had to be modified to model propulsion and fuselage effects on the aircraft stability. 
The stability derivatives used to validate the AVL model were from Ref. [34] in the cruise 
condition and are summarized in Table 6.1. The derivatives are listed in order of importance 
1 http: //people. nas.nasa.gov/ aftosmis/cart3d/cart3Dhome.html 
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to the aircraft dynamics [79]. The strategy used for matching results placed the greatest 
emphasis on matching stability derivatives from top to bottom. 

Table 6.1: Cessna 182T flight data stability derivatives at the cruise flight condition, angle 
units in radians [15]. 


Longitudinal 

Value 

Lateral/Directional 

Value 

c La 

4.41 

c h 

-0.0923 

C nla 

-0.613 

Cn 0 

0.0587 

Cm q 

-12.4 

Cl p 

-0.484 

c Lq 

3.9 

Cn p 

-0.0278 



c lr 

0.0798 


For comparison, the stability derivatives of a simple Cessna 182T geometry are presented 
first to give an indication of why the fuselage and propulsion effects needed to be modeled. 
This simple AVL model, consisting of only the lifting surfaces, is shown in Fig. 6.4. For 



the comparison of the simple AVL model to the flight data, only the longitudinal stability 
derivatives are shown in Table 6.2. 


Table 6.2: Comparison of Cessna 182T longitudinal stability derivatives between the simple 
AVL model and the flight data, angle units in radians [34], 



c La 

r 

^rricc 

r 


c LSe 

r 

K - y Tns e 

Flight Data 

4.41 

-0.613 

-12.4 

3.9 

0.43 

-1.122 

Simple AVL Model 

5.25 

-1.501 

-15.7 

9.5 

0.62 

-1.817 

Error 

19% 

144% 

27% 

143% 

44% 

70% 
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The large error in the pitch stiffness, C ma , and pitching damping, C mq , results in a 
conservatively stable aircraft where the pitch dynamics will be overly damped. Evidence of 
this can be seen in Ref. [80] , where the Cessna 182 dynamic model used in the perturbation 
and discrete gust responses was very stable. With very little control effort, the closed-loop 
system was able to alleviate the longitudinal disturbances that were described in Section 4.3. 
The static margin is calculated as the ratio of the pitch stiffness derivative over the lift curve 
slope. The large error in the pitch stiffness derivative indicates an overly large estimate of 
the static margin, and thus the location of the neutral point. The static margin of the Cessna 
182T used to obtain the stability derivative data was 14%, or 9.53 feet measured from the 
front of the propeller spinner. The AVL results for the simple geometry predicts a static 
margin of 29%, more than double that of the flight data. 

To obtain performance benefits by the reduction in static margin, a reasonable estimate 
of static margin from the aerodynamic tool was essential. Modeling the lifting surfaces alone 
clearly was insufficient in capturing the stability characteristics of the system. Discussed 
in Section 5.1, the fuselage and propulsion effects must be modeled in AVL to accurately 
represent the destabilizing effects they add. In Section 7.3 of Ref. [29], Etkin discusses the 
influence of the propulsion system on aircraft pitch stiffness. Running propeller alone has 
a destabilizing effect by reducing the pitch stiffness derivative through a radial force in the 
plane of the propeller. With the propeller in a tractor configuration, the slipstream creates 
a lift increment that is linear with alpha, effectively increasing the lift coefficient. This 
reduces the effect of the horizontal stabilizer on the neutral point location, decreasing the 
pitch stiffness and static margin. This effect will be even greater when using more than 
one propeller located directly in front of the wing. With the horizontal directly behind the 
propeller for the Cessna 182T, the downwash on the tail is altered and the effective angle of 
attack can be drastically reduced in the propeller wake. An example given by Etkin indicated 
a forward movement of the neutral point of 28% c, where c is the reference length [81]. 
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Since AVL was unable to directly model a propulsion system or its effects, the effects 
had to be simulated in the geometry definition. A canard-type lifting surface was added 
ahead of the horizontal fuselage surface to simulate the destabilizing effects of the propeller 
as seen in Fig. 5.6. Additionally, the section lift curve slopes on the main wing and the 
horizontal tail had to be adjusted. A default in AVL is to assume that each section lift curve 
slope is 27 r. Obviously, for a three-dimensional wing this is unrealistic and so the section 
lift curve slopes were adjusted to the 2-D airfoil experimental results from [82], However, 
this still produced an overly stable estimate of the neutral point. The size and location of 
the simulated propeller and the section lift curve slopes were used to adjust the AVL results 
to match the flight data stability derivatives. Since the horizontal stabilizer was directly 
behind the propeller, the section lift curve slopes were more heavily adjusted to simulate the 
reduced local angle of attack at that location. Finally, the gains on the elevator were adjusted 
to best match the control power derivatives for lift and pitching moment. Improving the 
accuracy of Cl s reduced the accuracy of C mg and so the elevator gain was chosen to split 
the two. Table 6.3 summarizes the comparison of the longitudinal derivatives. The flight 

Table 6.3: Comparison of Cessna 182T longitudinal stability derivatives from AVL to flight 
data, angle units in radians. 



c La 

c 

^ met 

r 

^rriq 

c Lq 

c L5e 

r 

Flight Data 

4.41 

-0.613 

-12.4 

3.9 

0.43 

-1.122 

Cessna 182T Model 

4.40 

-0.617 

-12.4 

5.9 

0.39 

-1.231 

Error 

-0.2% 

0.7% 

0% 

51% 

-9% 

8% 


data neutral point was 9.53 feet as measured from the front of the propeller spinner. With 
the adjustments described here, the predicted neutral point in AVL was 9.54 feet, a drastic 
improvement from the simple Cessna 182T AVL geometry. 

The pitch- rate dependent lift stability derivative, Cl , was unable to be matched. How- 
ever, according to Ref. [79] it is the least important longitudinal stability derivative, even 
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including stability derivatives not included in this analysis, such as the tuck derivative. Be- 
cause of the lack of influence on the longitudinal stability this error was considered to be 
irrelevant. 

Validation of the later al/directional stability derivatives was straightforward compared 
to the longitudinal stability derivatives. The vertical cross section of the fuselage was added, 
but the effects were over-destabilizing resulting in negative yaw stiffness (C nf3 < 0). The 
cross sectional area encompassing the windshield, side windows, and rear window was re- 
moved, improving the yaw stiffness. The influence of the vertical stabilizer on the dihedral 
derivative, Cip , and the yaw stiffness, C Ufj , was under-predicted with the bottom of the ver- 
tical ending at the top of the rear fuselage. This under-estimates the theoretical reference 
area and therefore the vertical stabilizer was extended down to increase the area, improv- 
ing the yaw stiffness and dihedral derivative. A summary of the AVL lateral/directional 
stability derivatives compared to the flight data values is given in Table 6.4. AVL under- 

Tablc 6.4: Comparison of lateral/directional stability derivatives from AVL to flight data, 
angle units in radians [34], 



c h 

Cn p 

Ci p 

c nr 


c lr 

Cy, 

Cy r 

C Y 

Flight Data 

-0.092 

0.059 

-0.484 

-0.094 

-0.028 

0.080 

-0.393 

0.214 

-0.075 

AVL Model 

-0.071 

0.057 

-0.449 

-0.107 

-0.020 

0.113 

-0.265 

0.259 

-0.046 

Error 

-23% 

-3% 

-7% 

-14% 

29% 

41% 

-33% 

21% 

-39% 


predicts the dihedral derivative by over 20% and unfortunately this is the most important 
lateral/direction stability derivative. The dihedral derivative is a function mainly of the 
geometry, specifically the wing dihedral, wing location (low, mid, high), and wing leading 
edge sweep angle. The Cessna 182T geometry has very low wing dihedral as discussed in 
Section 3.2. In a high-wing configuration at a non-zero sideslip angle, vortices form in the 
flow under with wing increasing the circulation and overall lift of the one wing [15]. This 
effect is stabilizing but is due to the viscous effects of the flow which cannot be captured in 
a potential flow vortex lattice code. This indicates that, clue to the high-wing Cessna 182T, 
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Table 6.5: Control derivatives from Cessna 182T AVL model and flight data, units in radi- 
ans [30]. 




Ci Sr 

c Ysr 

Cn Sa 

Cn Sr 

Flight Data 

-0.229 

0.0147 

0.187 

0.0216 

-0.0645 

Cessna 182T 

-0.231 

0.0155 

0.145 

0.00372 

-0.0688 

Error 

1% 

5% 

-22% 

-82% 

7% 


AVL will under predict the dihedral derivative. This was not corrected as it was desired to 
keep the AVL geometry as close to the actual geometry as possible. 

The lateral/directional control derivatives are summarized in Table 6.5. The gains for 
the aileron and rudder control surfaces were adjusted to best match the lateral/directional 
control derivatives. Three of the five control derivatives match well with one derivative 
having moderate error. The three most important control derivatives for a stability analysis 
match within 7% while C nSa completely misses the mark. Differences in drag between the 
asymmetrical deflection of the ailerons produces the yawing moment, but the drag due to 
the deflections was not accurately captured. The two derivatives that match poorly to the 
experimental data result as a secondary effect of the specific control surface deflection as 
opposed to the primary purpose of the control surface deflection. For example, a resulting 
side force is a secondary effect of a rudder deflection whereas a yawing moment is a primary 
effect. 

6.3.2 D8.2b Aerodynamic Model 

Experimental data for the D8.2b concept were unavailable for validation of the AVL 
model. Instead, Cart3D was used as a higher order method comparison [64,65]. A high- 
fidelity inviscid analysis software, Cart3D imports a triangulated surface mesh, generates 
a Cartesian volume mesh, executes a flow solution, and finally post-processes the results 
into forces and moments. A unique feature of Cart3D is the adjoint-based mesh refinement 
capability, which allows the user to converge a flow solution while minimizing error within 
the mesh. This eliminates the time consuming necessity of creating a mesh, performing 
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a mesh refinement study, and then modifying the mesh in areas of greatest error. A flow 
solution is generated for each mesh adaptation and was continued until the drag coefficient 
converged to within 5% of the previous adaptation cycle. 

The geometry used in Cart3D was taken from an OpenVSP model of the D8.2b, shown in 
Fig. 3.5, where a surface mesh was generated using the OpenVSP surface meshing capabilities 
exported as a triangulated file. The “CompGeom” functionality in OpenVSP computes the 
structured surface mesh and, given no open meshes are removed, guarantees a watertight 
geometry. The cell size of the structured mesh is governed by the number of interpolated 
cross sections in the model, which is controlled by the user. The structured mesh, exported 
as a triangulated surface mesh to Cart3D, is shown in Fig. 6.5. 



Figure 6.5: D8.2b mesh generated by CompGeom in OpenVSP. 

A sweep of angles of attack and sideslip at two Mach numbers was used to generate 
a dataset of aerodynamic data to compare with AVL; Mach 0.72 was used as the cruise 
condition and Mach 0.4 was used for the low-speed condition. A clean configuration — no 
slat or flap deflections with landing gear retracted — was used at both Mach numbers. All 
control surfaces were fixed in the neutral position. Low Mach numbers, i.e. Mach numbers 
where compressibility becomes insignificant, suffer from convergence issues in Cart3D, and 
therefore Mach numbers less than Mach 0.4 were not used. The angle of attack was varied 
from plus and minus six degrees in two degree increments at both Mach numbers. Sideslip 
angles were varied from plus or minus four degrees at a constant two degree angle of attack. 
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As mentioned in Section 5.1.2, the baseline D8.2b AVL model did not accurately capture 
the aerodynamics predicted by Cart3D. Modeling strategies employed on the Cessna 182T 
did not provide accurate results as would have been expected. As a result, an iterative process 
of modifying the baseline D8.2b AVL model was used to find the best modeling strategy to 
match the Cart3D results. Table 6.6 describes the iterative steps taken along with the 
abbreviations used in Figs. 6.6 and 6.7. Each model in the iterative process was the same 
as the previous with the exception of the single modification described. For example, only 
the horizontal fuselage cut airfoils were changed from the baseline model and all remaining 
components remained the same. 

Table 6.6: Description of D8.2b modeling steps with abbreviations defined. 
Abbreviation Description 

AVL Baseline AVL model closely matching OpenVSP model 

AVL Sym Horizontal fuselage cut airfoils made symmetric 
AVL Swing Horizontal fuselage removed, wing extended to centerline 
AVL Swing 0 Wing incidence angle removed 


Figure 6.6 shows the induced drag, lift, and pitching moment coefficient comparisons 
between Cart3D and the different AVL models of Table 6.6. For clarity, only Mach 0.72 
results are shown as the Mach 0.4 results follow the same trends. The induced drag coefficient 
of Fig. 6.6(a) shows that the different modeling strategies had very little effect on matching 
the lift-induced drag. At lift coefficients below approximately 0.8, AVL follows the trend of 
Cart3D until compressibility effects begin to dominate. Similar to the lift coefficient curve, 
an offset was present between Cart3D and AVL; due to the coupling between induced drag 
coefficient and lift coefficient, the offset in lift coefficient prediction results in the offset of the 
drag coefficient. Also, a zero lift pressure drag offset captured by Cart3D was not captured 
by AVL. The lift and pitching moment coefficients were very sensitive to the modeling 
strategy as indicated in Figs. 6.6(b) and 6.6(c). From Ref. [79], the lift and pitching moment 
slopes are the most important derivatives in the prediction of longitudinal stability; the 
greatest emphasis was placed on matching those derivatives to best predict the longitudinal 
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(a) Induced drag coefficient versus lift coefficient. 
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(b) Lift coefficient versus angle of attack. 
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(c) Pitching moment coefficient versus angle of attack. 


Figure 6.6: Drag, lift, and pitching moment coefficients for different AVL geometries. Data 
only from Mach 0.72 cases shown for clarity. 
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dynamics properties. The lift curve slope was impacted minimally by changes in the modeling 
geometry, but the lift coefficient at zero angle of attack changed greatly. AVL consistently 
over-predicted the zero angle of attack lift coefficient, with results improving as the horizontal 
fuselage was removed and the wing incidence removed. The lift curve slope was consistently 
under predicted when compared to Cart3D. Pitching moment coefficient was even more 
sensitive to the different models, especially when the horizontal fuselage was removed. With 
the horizontal fuselage, the pitching moment slope was drastically under predicted compared 
to Cart3D. The baseline model had a pitching moment offset at zero angle of attack and was 
the motivation for making the horizontal fuselage airfoils symmetric. With the symmetric 
airfoils for the fuselage, the zero angle of attack pitching moment offset was correct, but the 
under prediction in slope remained. Removing the horizontal fuselage reduced the pitching 
moment slope error, and removing the wing incidence decreased the offset at zero angle of 
attack. 

Lateral/directional force and moment curves are shown in Fig. 6.7. Removal of the 
horizontal fuselage improved the dihedral derivative, C) ;9 , while providing negligible impact 
on the yaw stiffness derivative, C np . Addition of the vertical fuselage resulted in an excellent 
agreement between Cart3D and AVL for the yaw stiffness, Fig. 6.7(b), while increasing error 
in the side force coefficient due to sideslip, C's 9 , seen in Fig. 6.7(c). The dihedral derivative 
was unaffected by the addition of the vertical fuselage. A compromise between the side 
force coefficient curve and the yaw stiffness had to be made as the addition of the fuselage 
improves one while degrading the other. As discussed by Stevens and Lewis in Ref. [79], 
the yaw stiffness derivative has greater impact on the aircraft dynamics than the side force 
coefficient. In fact, in Section 8.5 of Ref. [29], Etkin states, “...the whole derivative C\g\ p 
often has negligible effect on the vehicle dynamics.” Unfortunately, the dihedral derivative 
has some of the greatest impact on the flight dynamics and was unable to be matched with 
greater accuracy. As a function of dihedral angle, aspect ratio, and wing sweep [29], the 
dihedral derivative could not be improved as the geometry in terms of those parameters was 
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fixed. As such, the error had to be accepted as a limitation in the aerodynamic modeling 
capabilities of the vortex lattice. 

The baseline D8.2b AVL model did not accurately model the aerodynamics when com- 
pared to the Cart3D results. Figures 6.6 and 6.7 show that improvements can be made by 
modifying the AVL model. Using the comparison plots for the different D8.2b models, the 
final AVL model chosen for the multidisciplinary design optimization was when the horizon- 
tal fuselage lifting surface was removed, the simple wing extended to the centerline, wing 
incidence set to zero, and a vertical cut of the fuselage added. 

Force and moment coefficient comparison plots of the final D8.2b aerodynamic model at 
both Mach 0.72 and Mach 0.4 are shown in Fig. 6.8. With the horizontal fuselage removed 
and the vertical fuselage added, the agreement between Cart3D and AVL was quite good, 
especially in the longitudinal forces and moments. It is worth noting that the trailing edge of 
the fuselage contained a sharp edge allowing Cart3D to enforce a Kntta condition, resulting 
in the production of forces and moments from the fuselage alone. Were no sharp trailing 
edge present, the fuselage would produce no forces or moments, regardless of the wind 
vector orientation with respect to the fuselage body axis, as no Kntta condition could be 
applied. This is similar to performing an aerodynamic analysis of an ellipsoid in an Euler 
aerodynamic solver, where any angle of attack will result in zero net forces. This is due to 
the rear stagnation point having the freedom to align with the forward stagnation point, 
maintaining the same streamline as if the object was not in the fluid flow. 

The Mach number dependency is clearly visible for both Cart3D and AVL in all the 
plots of Fig. 6.8, with the exception of induced drag coefficient. At higher lift coefficients, 
compressibility drag begins to dominate as supercritical velocities begin to form on the 
surface, an effect captured by Cart3D. However, the compressibility correction in AVL, a 
Prandtl-Glauert correction [56], does not capture the increase in drag coefficient at Mach 
0.72 and lift coefficients greater than 0.85. At both Mach numbers, the induce drag coefficient 
was the same for AVL. In all the other force and moment curves of Fig. 6.8, the Mach number 
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Figure 6.7: Rolling moment, yawing moment, and side force coefficients for different AVL 
geometries. Data only from Mach 0.72 cases shown for clarity. 
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effects between Cart3D and AVL are similar, even though the slopes of the curves may differ. 
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Figure 6.8: Comparison of Cart3D and AVL force and moment coefficient predictions. 
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6.4 Cessna 182T Multidisciplinary Testing 


The Cessna 182T model was used to test the flight dynamics and system disturbance 
components to ensure each were behaving as expected. A simple trade study was performed 
by reducing the horizontal and vertical tail volume coefficients, described in Section 6.3.1. 
Results of the trade study are presented with the variation in drag coefficient versus the 
volume coefficient shown. Additionally, dynamic response plots and continuous turbulence 
power spectral density plots are given. 

Figure 6.9 shows the total drag coefficient as a function of the horizontal and vertical 
tail volume coefficients. As expected, the total drag decreased as the tail volume coefficients 
were reduced; the total drag was sensitive to both volume coefficients, but had a greater 
sensitivity to the vertical tail volume coefficient. Figure 6.10(a) shows the reduction in 
drag that resulted from decreasing the static margin. The drag benefit seen in Figs. 6.9 
and 6.10(a) comes from the reduction in wetted area and the resulting decrease in parasitic 
drag. When trimmed, induced drag for this geometry actually increased with the reduction 
of static margin as shown in Fig. 6.10(b). This explains why the total drag was more sensitive 
to varying the vertical tail volume coefficient as mentioned previously. The resulting increase 
in induced drag comes not from the reduction of static margin but the reduction in the span, 
and increase in lift coefficient, of the horizontal stabilizer. 

The elevator deflection angle, shown in Fig. 6.11, went from a slight negative deflection 
to a positive deflection. The Cessna 182T baseline geometry required less then a half degree 
of negative elevator deflection to trim at the cruise condition, which means the downwash 
from the main wing provided nearly sufficient induced angle on the horizontal stabilizer to 
trim. As the horizontal tail area and span were decreased, the required trim deflection on 
the stabilizer decreased, became neutral, and eventually positive. Even as the trim deflection 
increased the induced drag remained nearly constant, increasing by less than one percent. 
Because the location of the wing and center of gravity were fixed, the horizontal tail has to 
generate the same amount of moment to trim as the baseline case. As a result, a net increase 
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Figure 6.9: Cessna 182T drag coefficient for varying tail volume coefficients. 

in configuration induced drag was realized from the decreased horizontal stabilizer decreased 
efficiency. This indicates that simply reducing the static margin does not decrease induced 
drag but that a system level design, with static margin being allowed to decrease, must be 
performed in order to achieve the induced drag benefits of relaxed static stability. 

Initially, using moderate atmospheric disturbances with a probability of exeedance level 
of 1CT 3 , none of the configurations passed the heading hold turbulence check specified in 
SAE-AS94900. This resulted from the active control system not sufficiently driving the 
rudder in continuous turbulence. A great feature of the LQR is the flexibility in allowing 
for adjustments in the weighting matrices if a desired performance is not achieved. Due 
to the consistent failures in the heading hold check, the -0 weighting term in Eq. 4.29 was 
increased by a factor of 10, placing heavier penalties on non-zero values of yaw. This simple 
adjustment enabled all configurations to achieve adequate dynamic performance in the cruise 
flight condition. It should be noted, again, that the cruise flight condition alone does not 
adequately stress the control system for sizing the stabilizer surfaces. Additional flight 
conditions must be added to ensure proper sizing of the stabilizer surfaces. 
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(a) Drag coefficient vs. static margin. (b) Induced drag coefficient vs. static margin 

Figure 6.10: Sensitivity of total and induced drag coefficients to static margin. 



Figure 6.11: Elevator deflection angle versus static margin. 






Figures 6.12-6.17 show the baseline configuration’s dynamic response with the heavier 
yaw weighting to perturbations, continuous gust fields, and discrete gusts. These distur- 
bances were used to check the dynamic performance of the active control system as sum- 
marized in Table 4.3. The airspeed hold and longitudinal gust responses were omitted for 
brevity due to their similarity to the pitch hold and vertical gust responses. 

Time histories of the transient response to the perturbations specified by SAE-AS94900 
are shown in Figs. 6.12 and 6.13. The figures were separated into longitudinal and lateral 
states, showing both the short time response — less than ten seconds — and the long time 
response — total response over 100 seconds. All control surface states shown in Figs. 6.12, 
6.13, 6.16, and 6.17 are the total deflection angles, the sum of the perturbation and the 
steady-state deflection angle. Visible in the open- loop response of Fig. 6.12 is the phugoid 
mode which has been attenuated with little control power. Figure 6.13 shows the unstable 
spiral mode in the open-loop response that has been stabilized through a combination of 
aileron and rudder control surface deflections. 

Power spectral densities of the mean-squared responses are plotted in Figs. 6.14 and 6.15. 
With a natural frequency of 0.19 rad/s, the phugoid sensitivity can be seen in the vertical 
continuous turbulence responses as a spike in the mean-squared response of the open-loop 
system. Less obvious is the short period mode, which had an open-loop natural frequency of 
7.66 rad/s, that can be seen as a slight hump in the angle of attack, pitch rate, and pitch angle. 
The spikes in the lateral state responses of Fig. 6.15 correspond to the dutch roll natural 
frequency of 5.14 rad/s. The active control system successfully reduced the magnitude of the 
response of the phugoid and dutch roll modes while the short period response in turbulence 
was unchanged. 

Figures 6.16 and 6.17 show the response to a lateral and vertical discrete gust. As with 
the pitch perturbations, the phugoid is quite visible in the open-loop response. Due to the 
model being extremely stable at cruise, even in a discrete gust, the required elevator deflection 
of the closed-loop system was extremely small while providing a heavily damped response. 
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The closed-loop lateral gust response, Fig. 6.16, provides increased damping when compared 
to the open-loop response, as would be expected. Heavier weighting in the performance index 
on the heading angle, ip, improves the heading angle response but results in a sacrifice in the 
roll angle response. The heading perturbation state of Fig. 6.16 deviates from steady-state 
in response to the gust but immediately returns to the steady-state value, whereas the roll 
perturbation state has two full oscillation cycles before returning to the steady state. 
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Figure 6.12: Pitch hold perturbation check. 



























Figure 6.17: Vertical discrete gust response. 






Chapter 7 

Results and Discussion 


Five optimization cases, each with different constraints and design variables, were run 
in an exploration of the design space, searching for the active design constraints and sys- 
tem design sensitivities. All the optimization cases began using the design variables from the 
baseline configuration, Fig. 7.1, giving identical starting points for all optimized designs. The 
shorthand labels for the optimization cases with the applied constraints indicated are shown 
in Table 7.1, ranging from a nearly unconstrained optimization with traditional volume co- 
efficients to both static trim and dynamic response constraints applied to the optimization. 
The fixed tail volume coefficient optimization case used horizontal and vertical tail volume 
coefficients of 1.00 and 0.09 respectively, typical values for a jet transport [1]. Additionally, 
design cases with failed takeoff or minimum rate of climb capability at cruise were credited 
with zero range and fuel, resulting in a failed design as mission range was not met. The base- 
line configuration of Fig. 7.1 was unable to meet second segment climb gradient requirements 
as designed. 

Table 7.1: Shorthand labels of each optimization case with the applied constraints indicated. 


Case Label 

Fixed Tail 
Volume 

TOL 

Constraints 

Static 

Constraints 

Dynamic 

Constraints 

Fixed Static 
Margin 

FixedTailVol 

/ 

/ 



/ 

StaticConFixSM 


/ 

/ 


/ 

StaticConFreeSM 


/ 

/ 



SDynConF ixSM 


/ 

/ 

/ 

/ 

SDynConFreeSM 


/ 

/ 

/ 



Figure 7.1 shows a top and side view of the baseline configuration AVL model resulting 
from the aerodynamic validation of Section 6.3.2. As designed, the baseline configuration 
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was unable to meet the specified takeoff and landing constraints, with field length distances 
of 11,076 and 7,081 feet respectively. Federal Aviation Administration requirements state 
that a positive climb gradient must be achieved in the event of an engine failure during 
second segment climb, also failed by the baseline configuration. Tables 7.2 and 7.3 show 
the complete analyses of the static trim and dynamic response results with exceeded limits 
indicated in red. The constraints were unused in the analysis (the baseline configuration was 
fixed geometry and not optimized) and are shown for comparison to the cases with static 
trim and/or dynamic response constraints. The elevator effective angle of attack exceeds the 
constraints in both longitudinal perturbations in the stall condition. A residual roll angle 
greater than the allowed plus or minus one degree after five seconds was experienced in the 
cruise condition. 


Table 7.2: Baseline configuration static trim constraints. 


Condition 

Control Surface 

Deflection Angle 

OEI 

Aileron 

-0.33 


Rudder 

-19.64 

Takeoff 

Elevator 

1.65 

Maneuver 

Elevator 

8.32 


Figures 7.2 and 7.3 show the baseline configuration time response for the exceeded 
constraint limits highlighted in red in Table 7.3. The roll residual time responses in the 
cruise and stall flight conditions, Fig. 7.2, appear very similar. The roll angle starts at the 
perturbed value of five degrees and then decreases as the system returns to the steady-state 
condition. The horizontal dashed lines indicate the constraint limits and the vertical dashed 
line indicates the time that the roll angle must be within, and remain within, the constraint 
limits. The baseline configuration roll angle response in the cruise condition was too slow, 
failing to return to the allowed limits five seconds after the disturbance. In stall, the roll 
residual angle had an acceptable response. 
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Tabic 7.3: Baseline configuration dynamic response performance, in degrees unless otherwise 
noted. 


Condition 

Response Type 

Cruise Condition 

Stall Condition 

Discrete Gust 

Aileron 

0.00 

5.88 


Rudder 

-1.38 

-1.88 


Elevator (Long) 

3.17 

8.42 


Elevator (Vert) 

3.23 

9.52 

Turbulence 

RMS Pitch 

0.05 

1.08 


RMS Roll 

1.29 

1.16 


RMS Yaw 

0.36 

0.42 

Perturbation 

Aileron 

0.00 

6.37 


Rudder 

4.46 

1.75 


Elevator (Pitch) 

7.41 

20.49 


Elevator (Speed) 

5.45 

24.37 

Residual 

Pitch 

-0.24 

-0.29 


Roll 

-1.77 

-0.64 


Speed (fps) 

0.00 

0.00 


The longitudinal perturbation responses of the baseline configuration are shown in 
Fig. 7.3. In the stall condition, the elevator was heavily loaded, almost to the upper con- 
straint, prior to any perturbations being introduced into the system. When perturbed in 
pitch and airspeed, the elevator responded to the system disturbances and exceeded the al- 
lowable elevator effective angle of attack. This high loading during the perturbation response 
would cause the horizontal stabilizer to enter the nonlinear aerodynamics region, potentially 
entering a full stall and risking the loss of control. 

Table 7.4 compares the design variables from each optimization to the baseline configu- 
ration. Side and top view comparisons of the AVL geometries are given in Figs. 7.4 and 7.5. 
Between all the optimal designs, thrust and wing area varied minimally as these were driven 
by the takeoff and landing held constraints. The only exception was the SDynConFixSM 
case with the wing area increasing to over 1,700 square feet. The wing apex location, for 
all designs, shifted forward, increasing the horizontal stabilizer moment arm. This was most 
evident in the optimization case where the traditional tail volume coefficient was used. Re- 
ducing horizontal tail area resulted in an empty weight and wetted surface area reduction, 
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(a) Top view. 


(b) Side view. 


Figure 7.1: Top and side view of baseline D8.2b configuration, modeled in AVL. 

a property exploited by the optimizer in fixed tail volume coefficient configuration. As the 
stabilizers’ areas were calculated based upon fixed volume coefficients, the moment arms 
for both the vertical and horizontal stabilizers were maximized by shifting the wing apex 
to the forward- most limit and the stabilizers to the rear- most limit. This case highlights 
the limitation of multidisciplinary design optimization when only traditional performance 
requirements are used, and stability and control considerations are only captured through 
the inclusion of a volume coefficient. 

Vertical stabilizer area varied in size throughout the five optimization cases. The ver- 
tical stabilizer area reached a minimum in the StaticConFreeSM case where the design was 
given the greatest flexibility without the burden of the dynamic response constraints. When 
the dynamic response constraints were added, the vertical stabilizer area increased to near 
the baseline configuration in the SDynConFixSM case. Having static margin as a design 
variable allowed the vertical stabilizer area to decrease, considering it a weight and vis- 
cous drag penalty. A similar behavior was seen with the horizontal stabilizer. Again, in 
the StaticConFixSM case, the horizontal stabilizer area increased to slightly more than the 
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Figure 7.2: Baseline configuration cruise and stall roll residual time responses. Horizontal 
dashed lines indicate the roll residual constraint limits and the vertical dashed line indicates 
the time where the residual constraint became active. 

baseline configuration, but freeing static margin allowed the optimizer to freely shift the 
center of gravity and neutral point, reducing the required horizontal stabilizer area. 

Static margin was allowed to vary for two designs, and in both cases the optimizer 
chose to increase the static margin, opposite as what would have been expected as discussed 
previously in Section 2.1. As a dependent variable, the center of gravity was placed based 
upon the static margin, a design variable. The optimizer discovered, as discussed in the 
detailed description of each configuration, a benefit of shifting the center of gravity forward, 
increasing the tail moment arm and providing greater control authority to the stabilizers. 
This allowed the stabilizers to be decreased in size giving a weight and viscous drag benefit. 
This was a unique solution in a design space not typically explored with a fixed center of 
gravity position; increasing the open loop stability provided a greater reduction in system 
fuel burn, through reduced stabilizer sizing, than the benefits of relaxed static stability for 
the given design space and controller design. 
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(a) Pitch perturbation elevator response. (b) Airspeed perturbation elevator response. 


Figure 7.3: Baseline configuration elevator time responses to pitch and airspeed perturbations 
in the stall condition. The horizontal dashed lines indicate the elevator effective angle of 
attack limits. 



(a) Baseline configuration. 



(c) StaticConFixSM configuration. 


(e) SDynConFixSM configuration. 



(b) FixedTailVol configuration. 



(d) StaticConFreeSM configuration. 



(f) SDynConFreeSM configuration. 


Figure 7.4: Side view of baseline configuration and all optimization cases. 
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(a) Baseline configuration. (b) FixedTailVol configura- (c) StaticConFixSM configuration. 

tion. 



(d) StaticConFreeSM configura- (e) SDynConFixSM configuration, (f) SDynConFreeSM configuration, 
tion. 


Figure 7.5: Top view of baseline configuration and all optimization cases. 
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7.1 Optimal Design using Fixed Tail Volume Coefficients 


The optimal design using fixed tail volume coefficients with takeoff and landing (TOL) 
field constraints is shown in Fig. 7.6, and the associated design variables given in Table 7.5. 
As mentioned previously, the wing apex position was moved to the forward limit, 20 feet 
measured from the fuselage nose, while the vertical stabilizer sweep was increased to the 
upper limit, 65 degrees, shifting the horizontal tail to the aft limit. For a fixed horizontal 
tail volume coefficient, increasing the tail moment arm results in a reduction in horizontal 
stabilizer area, and thus wetted area, ergo viscous drag. With the vertical stabilizer root 
chord fixed to the fuselage rear and sweep limited to 65 degrees, the moment arm between the 
wing and vertical stabilizer aerodynamic centers was constrained, and for fixed tail volume 
coefficient caused the vertical stabilizer area to increase. Wing sweep was increased, giving 




(b) Side view. 


Figure 7.6: FixedTailVol configuration top and side views. 


Table 7.5: FixedTailVol design variable summary. 


T 

S 

A 

T 

S H t 

SyT 

A VT 

X w 

»' apex 

SM 

(lb) 

(ft 2 ) 

(deg) 

(deg) 

(ft 2 ) 

(ft 2 ) 

(deg) 

(ft) 

(%) 

22,590 

1,451 

29.5 

0.1 

222.1 

113.3 

65.0 

20.0 

10 


a structural weight penalty, to take advantage of compressibility drag reductions at the cruise 
condition, resulting in fuel burn savings. Additionally, the wing sweep shifted the neutral 
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point aft, and thus the center of gravity aft, increasing the nose up pitching moment due to 
the wing aerodynamic center being forward of the CG. 

Figure 7.7 shows the takeoff and landing field performance of the design optimization 
using specified tail volume coefficients. Takeoff and landing were the only constraints on 
the optimization, other than mission range had to be met with a successful takeoff. Using 
only the best design from each generation, Fig. 7.7(a) shows the takeoff and landing field 
performance for each design, normalized by the constraint field lengths of 8,000 and 5,800 
feet, respectively. To determine how active a constraint was throughout the design optimiza- 
tion, the takeoff and field length performance for every design evaluated in the optimization 
is shown in Fig. 7.7(b) as a histogram. Indicated in the histograms, takeoff and landing 
field length were active constraints in the design space with normal distributions centered 
near each field length constraint. Wing area and engine static thrust sizes are dependent on 
takeoff and landing field requirements, and any decrease in held length results in increased 
weight due to increased wing area and/or engine static thrust. No fuel burn benefit results 
from held length performances better than the constraints, so the optimizer minimized wing 
area and static thrust, pushing the design to the constraints. 
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(a) Normalized takeoff and landing field lengths. 


(b) Histogram of takeoff and landing field distances. 


Figure 7.7: FixedTailVol configuration takeoff and landing held length performance shown 
normalized by the constraints (a) and held length distances of all evaluated designs in a 
histogram (b). 
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The static trim and dynamic response analyses were run for the FixedTailVol case with 
results given in Tables 7.6 and 7.7. These constraints were not applied to the FixedTailVol 
optimization case, only takeoff and landing held constraints were applied, but the results 
give an indication of the shortcomings of only using performance constraints in an opti- 
mization. Red font in Tables 7.6 and 7.7 indicates a value that would have exceeded an 
allowable constraint limit had it been applied. In the one engine inoperative condition, the 
rudder deflection required to trim the yawing moment would result in stalling of the surface. 
The stall condition dynamic response analysis results of Table 7.7 show several limits that 
were exceeded. If the static trim and dynamic response constraints had been applied the 
FixedTailVol design would have been declared an infeasible design, rejecting it from the de- 
sign space. Ignoring stability and control during the optimization, i.e., sizing the stabilizers 
using only tail volume coefficients, resulted in a design incapable of meeting all the handling 
and controllability requirements. 

Table 7.6: FixedTailVol configuration static trim defiections, degrees. 


Condition 

Response Type 

Deflection 

OEI 

Aileron 

0.19 


Rudder 

-33.26 

Takeoff 

Elevator 

0.86 

Maneuver 

Elevator 

10.89 


The elevator effective angle of attack time responses to pitch and airspeed perturbations 
are shown in Fig. 7.8. After the perturbation was introduced into the system at time zero, 
the elevator effective angle of attack increased to approximately 20 degrees, far exceeding 
the allowable loading on the horizontal stabilizer. As the system returns to the steady-state 
condition after four seconds, the trim effective angle of attack was near the allowable limit, 
which indicated that nearly the maximum allowed loading was required to trim in the stall 
flight condition. With the stabilizer fully loaded for trim in the stall condition, no additional 
margin was available to respond to a system disturbance. 
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Tabic 7.7: FixedTailVol configuration dynamic response performance, in degrees unless oth- 
erwise noted. 


Condition 

Response Type 

Cruise Condition 

Stall Condition 

Discrete Gust 

Aileron 

0.00 

-11.79 


Rudder 

1.92 

7.40 


Elevator (Long) 

1.00 

11.06 


Elevator (Vert) 

1.12 

11.33 

Turbulence 

RMS Pitch 

0.06 

0.52 


RMS Roll 

0.49 

1.64 


RMS Yaw 

0.40 

2.63 

Perturbation 

Aileron 

0.00 

8.41 


Rudder 

7.80 

5.88 


Elevator (Pitch) 

5.51 

18.28 


Elevator (Speed) 

3.42 

22.33 

Residual 

Pitch 

-0.18 

-0.42 


Roll 

-1.43 

-0.61 


Speed (fps) 

0.00 

0.00 


Figure 7.9 shows the FixedTailVol configuration roll response to a five-degree system 
perturbation in the cruise and stall flight conditions. After the disturbance, the system begins 
to return to the steady-state condition, but overshoots, and has to return. The overshoot in 
the cruise flight condition slowly returns to zero but fails to achieve the allowable residual 
five seconds after the disturbance. The stall condition roll residual behavior exhibited a 
slightly different behavior. The initial overshoot hits the lower bound of the residual limit, 
but returns to within the allowed residual prior to the five second cutoff. Unlike the cruise 
condition, a smaller, second oscillation occurred but never exceeded the limits. 
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(a) Pitch perturbation elevator response. (b) Airspeed perturbation elevator response. 

Figure 7.8: FixedTailVol configuration elevator time response to pitch and airspeed per- 
turbations in the stall flight condition. The horizontal dashed lines indicate the maximum 
allowable elevator effective angles of attack. 




(a) Cruise roll residual. (b) Stall roll residual. 

Figure 7.9: FixedTailVol configuration roll residual time responses in the cruise and stall 
flight conditions. The horizontal dashed lines indicate the maximum allowed deviations 
after five seconds, and the vertical dashed line indicates when the residual response was 
applied. 
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7.2 Optimal Designs with Static Trim Constraints 


Applying static trim constraints to the multidisciplinary design optimization drastically 
altered the optimal designs compared to the case using traditional volume coefficients. Fig- 
ures 7.4 and 7.5 compare the optimal designs with static trim constraints, both fixed and 
free static margin, with all the optimal design cases. Clearly, the vertical tail area and sweep 
were reduced from the FixedTailVol case, and the wing apex location, especially in the fixed 
static margin case, was moved aft toward the center of the fuselage. Optimizing with the 
static trim constraints instead of fixed tail volume coefficient resulted in very little change 
in wing area, engine static thrust, and wing sweep as indicated in Table 7.4, but a nonzero 
dihedral angle was introduced. As the tail areas were sized using static trim constraints, the 
horizontal and vertical stabilizer areas decreased. This indicated using fixed tail volume coef- 
ficients resulted in a poor optimal solution because the stabilizer areas were greater than the 
cases with static trim constraints, and the one engine inoperative trim deflection resulted in 
a stalled control surface. Removing the fixed tail volume coefficient requirement and adding 
static trim constraints allowed the wing apex to shift aft, allowing for a configuration with 
reduced load requirements on the stabilizers. 

Allowing static margin to be a free design variable resulted in several design changes 
as seen in Fig. 7.10. The two configurations have a similar wing planform with the 
StaticConFreeSM wing apex farther forward than the StaticConFixSM case. Static mar- 
gin increased to 15.1%, allowing for reduced horizontal and vertical stabilizer areas and 
wetted area reductions. Having only 10.6 degrees of vertical tail sweep results in structural 
weight and wetted area reductions. Table 7.8 summarizes the design variables of the two 
configurations. 

Similar to Section 7.1, takeoff and landing field lengths were active constraints in the 
design space, heavily driving the wing sizing and static thrust. As gross weight decreased 
from the FixedTailVol configuration, both wing area and engine static thrust could be de- 
creased to meet the same field length requirements. Figure 7.11 shows the normalized takeoff 
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(a) StaticConFixSM top view. 


(b) StaticConFixSM side view. 



(c) StaticConFreeSM top view. 



(d) StaticConFreeSM side view. 


Figure 7.10: StaticConFixSM and StaticConFreeSM optimization cases top and side views. 


and landing constraints along with the associated histograms for all designs analyzed during 
the genetic algorithm search. The optimization pushed the upper bounds of the constraints 
as the histogram shows a large number of designs exceeding the 8,000 and 5,800 feet held 
lengths. In the case of takeoff, numerous designs had held lengths less than 8,000 feet, but 
at a greater system total fuel burn, increasing the system design objective function. 


Table 7.8: Static trim cases design variable summary. 


Design 

T 

S 

A 

r 

S H t 

Svt 

A VT 

X w 

v* apex 

SM 


(lb) 

(ft 2 ) 

(deg) 

(deg) 

(ft 2 ) 

(ft 2 ) 

(deg) 

(ft) 

(%) 

StaticConFixSM 

22,200 

1,422 

30.0 

8.8 

118.8 

50.4 

26.9 

35.6 

10.0 

StaticConFreeSM 

22,120 

1,411 

28.5 

8.1 

102.1 

32.4 

10.6 

23.4 

15.1 
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Distance, ft 


(a) StaticConFixSM normalized takeoff and landing (b) StaticConFixSM histogram of takeoff and landing 
field lengths. field lengths. 
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(c) StaticConFreeSM configuration normalized take- (d) StaticConFreeSM configuration histogram of 
off and landing field lengths. takeoff and landing field lengths 


Figure 7.11: StaticConFixSM and StaticConFreeSM configurations takeoff and landing field 
length performance shown normalized by the constraints and field length distances of all 
evaluated designs in a histogram. 
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The static trim constraints included one engine inoperative (OEI) aileron and rudder 
deflection, takeoff elevator effective angle of attack, and maneuver condition elevator effective 
angle of attack limits. As expected, OEI rudder deflection and takeoff elevator effective angle 
of attack, the effective angle of attack required to lift the nose wheel at rotation speed, were 
increasingly active as the optimization progressed, indicating the desire to minimize the 
horizontal and vertical stabilizer areas. Figures 7.12(a) and 7.12(b) show the normalized 
trim constraints of the best designs. Only the rudder and elevator constraints were active, 
especially in later generations, and aileron deflection angle was minimal during the one engine 
inoperative flight condition. 

Numerous designs had takeoff and maneuver trim constraints within the specified limits 
of -11.4 and 12.2 degrees. Figures 7.12(c) and 7.12(d) show the large distribution of elevator 
effective angles of attack. The maneuver elevator range of values were smaller in the fixed 
static margin case than the free static margin, but in both optimizations the majority of 
designs fell below the upper constraint of 12.2 degrees. Even though the maneuver elevator 
constraint was inactive in the best design of each generation, it can be seen that the constraint 
was active in the free static margin design space as several hundred designs violated the upper 
constraint bounds. Many designs evaluated in the fixed static margin case had maneuver 
elevator effective angles of attack approaching the upper limit, but few designs ever exceeded 
the upper limit. 

An interesting feature of the configurations was the positive elevator lift required to trim 
during the maneuver condition. The wing, as designed in Ref. [10], had airfoils shaped to 
reduce the nose down pitching moment. As the airfoils were again used in this research, the 
combination of a heavily loaded wing, the wing lift vector forward of the center of gravity, 
and a small wing aerodynamic pitching moment resulted in a net nose up pitching moment 
required to be counteracted by the elevator. As the geometric deflection angle of the elevator 
was negative, the combination of geometric angle of attack, downwash, and deflection angle 
resulted in a positive lift coefficient to counteract the nose up pitching moment of the heavily 
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Generation 



Generation 


(a) StaticConFixSM normalized deflection angles, (b) StaticConFreeSM configuration normalized de- 
degrees. flection angles, degrees. 
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(e) StaticConFixSM OEI deflection angles his- 
togram, degrees. 



Deflection Angle, deg 


(f) StaticConFreeSM configuration OEI deflection 
angles histogram, degrees. 


Figure 7.12: Static trim constraints with each generation best design deflection angles nor- 
malized by the constraints (a,b) and deflection angles histograms of evaluated designs (c-f). 
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loaded wing. Although the effective angle of attack may seem counterintuitive, geometrically 
the geometry was behaving as expected. This highlights the benefit of using the elevator 
effective angle of attack as the horizontal stabilizer sizing constraint instead of a geometry 
deflection angle. 

Static trim constraints were applied to the optimization and were used to size the 
horizontal and vertical stabilizers in parallel with the planform geometry. Including the 
constraints in the optimization eliminated the necessity of tail volume coefficients while still 
considering trim controllability. Table 7.9 shows that, unlike the FixedTailVol case, all the 
static trim constraints were met, with rudder and elevator deflection angles being the most 
active constraints. The stabilizer and control surfaces were adequately sized to provide trim 
without exceeding the allowed deflections constraints or stall the surface. 


Table 7.9: Trim deflection angles, in degrees, for static trim optimization cases. 


Flight Condition 

Control Surface 

StaticConFixSM 

StaticConFreeSM 

OEI 

Aileron 

0.45 

0.76 


Rudder 

-19.22 

-19.80 

Takeoff 

Elevator 

-11.25 

-11.38 

Maneuver 

Elevator 

4.26 

1.58 


Reviewing Fig. 7.10, the horizontal and vertical stabilizer areas, even in a twin vertical 
tail configuration, were very small. Using only static trim constraints does not adequately 
size the stabilizer surfaces as the dynamic response performance was far from acceptable. 
Table 7.10 highlights the poor dynamic performance of the elevator in the stall condition. 
The pitch perturbation maximum elevator deflection and residual pitch angle both exceed 
the constraints. Additionally, both the fixed and free static margin cases had deficient roll 
residual performance in the cruise condition, not accounted for during the optimization with 
static trim constraints only. 

Figures 7.13 to 7.15 show the time response plots for the constraints highlighted in 
red in Table 7.10. Figure 7.13 shows the elevator effective angle of attack time response 
to the airspeed perturbation for the StaticConFixSM and StaticConFreeSM configurations. 
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Tabic 7.10: Dynamic response performance for static trim optimization cases, in degrees 
unless noted otherwise. 


Condition 

Response Type 

StaticConFixSM 

Cruise 

Stall 

StaticConFreeSM 

Cruise 

Stall 

Discrete Gust 

Aileron 

0.00 

-11.55 

0.00 

-10.83 


Rudder 

-1.76 

8.44 

-1.62 

7.86 


Elevator (Long) 

-0.74 

4.89 

-1.70 

2.17 


Elevator (Vert) 

-0.88 

6.01 

-1.85 

3.54 

Turbulence 

RMS Pitch 

0.11 

0.74 

0.14 

0.81 


RMS Roll 

1.97 

2.33 

1.81 

2.05 


RMS Yaw 

0.37 

1.41 

0.37 

1.51 

Perturbation 

Aileron 

0.00 

6.72 

0.00 

6.89 


Rudder 

3.15 

3.77 

3.47 

3.80 


Elevator (Pitch) 

5.35 

13.60 

4.68 

11.16 


Elevator (Speed) 

2.56 

18.49 

-1.79 

16.21 

Residual 

Pitch 

-0.37 

-0.67 

-0.35 

-0.60 


Roll 

-1.14 

-0.77 

-1.34 

-0.75 


Speed (fps) 

0.00 

0.00 

0.00 

0.00 


An elevator effective angle of attack of far greater than the allowed 12.2 degrees occurred 
in both configurations at the stall flight condition, settling out six seconds after the initial 
airspeed perturbation. 

The residual pitch angle time responses for the StaticConFixSM and StaticConFreeSM 
configurations in the stall condition are shown in Fig. 7.14. Horizontal dashed lines indicate 
the maximum allowable pitch angle residual after five seconds, and the vertical dashed line 
indicates when the residual constraint began to be applied. The pitch angle response over- 
shoots the steady-state condition before returning, exceeding the allowable deviation after 
five seconds. At five seconds the StaticConFixSM had a pitch residual of -0.67 degrees, and 
the StaticConFreeSM had a slightly better response with a residual of -0.60 degrees. 

Table 7.10 shows that the static trim optimization cases failed to meet the roll residual 
constraint limits in the cruise flight condition. Figure 7.15 shows the roll angle residuals for 
both static trim constraint optimization cases in the cruise flight condition. Both cases failed 
to meet the allowed plus or minus one degree roll residual deviation from the steady-state 
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(a) StaticConFixSM configuration. (b) StaticConFreeSM configuration. 

Figure 7.13: StaticConFixSM and StaticConFreeSM configurations elevator response to an 
airspeed perturbation in the stall condition. The horizontal dashed lines indicate the maxi- 
mum allowed elevator effective angle of attack. 

condition five seconds after a perturbation. The responses between the two configurations 
were very similar, with the free static margin case having poorer performance. 
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(a) StaticConFixSM configuration. (b) StaticConFreeSM configuration. 

Figure 7.14: StaticConFixSM and StaticConFreeSM configurations pitch residual time re- 
sponse in the stall condition. The horizontal dashed lines indicate the maximum allowable 
residual pitch angle five seconds after the initial perturbation, and the vertical dashed line 
indicates when the residual response was applied. 




(a) StaticConFixSM configuration. (b) StaticConFreeSM configuration. 

Figure 7.15: StaticConFixSM and StaticConFreeSM configurations roll residual time re- 
sponses in the cruise condition. The horizontal dashed lines indicate the maximum allowed 
deviations after five seconds, and the vertical dashed line indicates when the residual response 
was applied. 
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7.3 Optimal Designs Using Static Trim and Dynamic Response Constraints 


Referring back to Figs. 7.4 and 7.5, comparisons can be drawn between the configura- 
tions optimized with and without dynamic constraints. The vertical stabilizer area increased 
for both the SDynConFixSM and SDynConFreeSM configurations, Figs. 7.4(e) and 7.4(f), 
compared to the StaticConFixSM and StaticConFreeSM configurations. Specifics of the dy- 
namic constraint optimized designs can be seen in Table 7.4 where the vertical stabilizer area 
was increased from the baseline for the fix static margin configuration. With static margin 
fixed, the total fuel burn was greatest of all the optimized designs, totaling 30,304 pounds. 
Freeing static margin and optimizing with static trim and dynamic response constraints re- 
duced the fuel burn by approximately 4,000 pounds. The SDynConFreeSM configuration 
gross weight, 156,767 pounds, was reduced to less than the baseline and FixedTailVol de- 
signs. Only the static constraint optimized designs, StaticConFixSM and StaticConFreeSM, 
had smaller gross and fuel weights. 

To compare the fixed versus free static margin designs directly, Fig. 7.16 shows the top 
and side views of the two configurations, with design variable details provided in Table 7.11. 
Wing area of the SDynConFixSM case was 1,730 square feet, the only optimized design 
that was not in the region of 1,450 square feet. Static thrust was less for the free static 
margin case while also having a wing area of 1,426 square feet. Figure 7.17 shows that the 
field length for the SDynConFreeSM case was active throughout the optimization, unlike the 
SDynConFixSM case. The larger wing area of the SDynConFixSM design, with larger thrust 
than the other designs, alleviated the field length constraint. The histogram of Fig. 7.17(b) 
shows that for all the designs evaluated during the optimization, rarely was the takeoff or 
landing constraint violated for the SDynConFixSM configuration. 

Takeoff and landing field lengths were active constraints for the SDynConFreeSM case, 
as indicated by the normalized field length plot and histogram of Figs. 7.17(c) and 7.17(d). 
Many of the evaluated designs had takeoff and landing field distances greater than the allowed 
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(a) SDynConFixSM top view. 


(b) SDynConFixSM side view. 


r 


■ 




(c) SDynConFreeSM top view. (d) SDynConFreeSM side view. 

Figure 7.16: SDynConFixSM and SDynConFreeSM optimization cases top and side view. 


limits, resulting in infeasible designs. The normalized takeoff and landing distances for the 
best designs in each generation consistently pushed the upper constraint bounds. 

Application of a moderate lateral discrete gust was used to stress the lateral/directional 
control system in the cruise and stall flight conditions. In the cruise flight condition for both 
the SDynConFixSM and SDynConFreeSM designs, deflection angles for both the aileron 
and rudder were small, less than two degrees for nearly all the designs, as indicated in the 


Table 7.11: SDynConFixSM and SDynConFreeSM configurations design variable summary. 


Design 

T 

S 

A 

r 

Sht 

SvT 

A VT 

A w 

vv apex 

SM 


(lb) 

(ft 2 ) 

(deg) 

(deg) 

(ft 2 ) 

(ft 2 ) 

(deg) 

(ft) 

(%) 

SDynConFixSM 

23,670 

1,730 

10.8 

9.7 

288.7 

79.2 

16.7 

51.5 

10.0 

SDynConFreeSM 

22,200 

1,426 

24.1 

9.4 

169.5 

44.5 

12.3 

43.2 

22.4 
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1.6 




Generation 


Distance, ft 


(a) SDynConFixSM normalized takeoff and landing (b) SDynConFixSM histogram of takeoff and landing 
field lengths. field lengths. 



0 10 20 30 40 50 60 70 80 90 100 110 

Generation 

(c) SDynConFreeSM normalized takeoff and landing 
field lengths. 



Distance, ft 


(d) SDynConFreeSM histogram of takeoff and land- 
ing field lengths. 


Figure 7.17: SDynConFixSM and SDynConFreeSM configurations each generation best de- 
sign normalized takeoff and landing field lengths and field length histograms of each evaluated 
design. 


histograms of Fig. 7.18. The stall flight condition was more demanding on the flight control 
system, but never becoming an active constraint in the best designs. The maximum aileron 
and rudder deflections were nearly constant for the SDynConFixSM best designs, whereas the 
maximum aileron and rudder deflections in the SDynConFreeSM case increased throughout 
the optimization, settling only when the best solution was found. The histograms show 
that the aileron response to a lateral gust consistently had the largest deflection in both 
optimization cases, but no designs ever exceeded the upper constraint limit of 20 degrees. 
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Figure 7.18: SDynConFixSM and SDynConFreeSM configurations aileron and rudder re- 
sponse to a lateral gust. 


The maximum allowed root-mean-squared (RMS) deviation from moderate continu- 
ous turbulence was 5 degrees in both pitch and yaw, and 10 degrees in roll. Neither the 
SDynConFixSM nor SDynConFreeSM optimization cases struggled with this constraint, 
Fig. 7.19, as the best designs never neared the constraints. The closest approach to a 
continuous turbulence constraint was yaw in the stall flight condition. The RMS yaw was 
the only continuous turbulence response that exceeded three degrees, occurring in very few 
designs. Pitch and roll RMS angles, at both stall and cruise flight conditions, were well 
below the constraints. 
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(c) SDynConFreeSM normalized RMS turbulence re- (d) SDynConFreeSM histogram of RMS turbulence 
sponse. response, degrees. 


Figure 7.19: SDynConFixSM and SDynConFreeSM configurations RMS turbulence response. 


Maximum aileron and rudder deflections angles from the five degree roll perturbation 
were small in the cruise condition, with many of the aileron deflections being near zero, visible 
in Fig. 7.20. As the histograms of Figs. 7.20(b) and 7.20(d) show, the cruise roll perturbation 
maximum rudder deflection for nearly all the designs evaluated was less than eight degrees, 
indicating that the roll perturbation in the cruise condition had very little influence on the 
design. With the small maximum aileron deflections, the closed-loop system was relying on 
the natural roll stability provided by the wing dihedral. 

In the stall flight condition, maximum aileron deflections were significantly higher than 
the cruise condition, with many designs having deflections between 12 and 14 degrees for 
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(a) SDynConFixSM normalized roll perturbation re- (b) SDynConFixSM histogram of roll perturbation 
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(c) SDynConFreeSM normalized roll perturbation re- (d) SDynConFreeSM histogram of roll perturbation 
sponses. responses. 


Figure 7.20: SDynConFixSM and SDynConFreeSM configurations roll perturbation response 
performance. 


both the SDynConFixSM and SDynConFreeSM cases. The rudder deflections in the stall 
condition closely resembled those in the cruise condition for both cases. Even though the 
aileron deflections were significantly higher in the stall condition, Figs. 7.20(a) and 7.20(c) 
show that the maximum deflections were far below the maximum of 20 degrees. Again, in 
terms of maximum deflections, the roll perturbation was inactive in the design space even 
for the stall flight condition. 

A second requirement on the system perturbation response was the time to return 
to steady-state conditions within a residual tolerance. Five seconds after an initial system 
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perturbation, the pitch and roll Euler angle residuals were required to return to within 0.5 and 
1.0 degrees of the steady state, respectively. The Baseline, FixedTailVol, StaticConFixSM, 
and StaticConFreeSM configurations all failed to have an adequate roll response with a roll 
residual greater than the allowed limit. Also, during preliminary optimizations with the 
static trim and dynamic response constraints, the roll residual constraint was consistently 
active, resulting in the genetic algorithm struggling to find feasible designs. Given this 
experience, the roll angle weighting of Eq. 4.36 was increased by a factor of 10. 

Figure 7.21 shows the normalized residual of the best designs and the residual histograms 
of all evaluated designs with the increased roll weighting. For both cases, neither the cruise 
nor the stall flight condition pitch residuals were constraining on the design space. Nearly 
all evaluated designs had pitch residuals within the allowed limits. In previous optimizations 
of the SDynConFixSM and SDynConFreeSM cases, a large number of designs exceeded the 
allowed roll residual in both the cruise and stall flight conditions. The normalized residuals 
of the best designs, shown in Figs. 7.21(a) and 7.21(c), show that the increased roll weighting 
reduced the roll residuals, removing it as an active constraint in the design space. Reviewing 
the histograms of Figs. 7.21(b) and 7.21(d), the heavier roll weighting made it to where there 
were no designs that exceeded the plus or minus one degree roll angle limit. 

The airspeed deviation after 30 seconds in all cases was driven to zero as seen in the 
normalized residual plots. If an altitude state with a fixed altitude requirement was placed on 
the linear simulation, this constraint could become active, but would also require a throttle 
model, not currently included in the system dynamics used in this research. As this residual 
constraint was always zero, it will be discussed no further. 

Elevator response to moderate longitudinal and vertical gusts is shown in Fig. 7.22. A 
large range of elevator effective angle of attacks resulted in the design space, none exceeding 
the constraint limits of -11.4 and 12.2 degrees. The best designs never approached the 
constraints while the center of the normal distributions in the histograms, the approximate 
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Figure 7.21: SDynConFixSM and SDynConFreeSM configurations perturbation residual re- 
sponses. 


average deflections, were greater in magnitude for the fixed static margin case in the stall 
flight condition. 

Contrary to the elevator discrete gust response, the maximum effective elevator de- 
flections in response to perturbations were large, with many designs exceeding the allowed 
effective angle of attack as indicated in Fig. 7.23. The airspeed perturbation was the most 
constraining disturbance in the sizing of the horizontal stabilizer, with many designs signif- 
icantly exceeding the upper limit of 12.2 degrees, shown in Figs. 7.23(b) and 7.23(d). The 
normalized effective angle plots indicate that the initial best designs violated this constraint, 
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Figure 7.22: SDynConFixSM and SDynConFreeSM configurations elevator responses to both 
longitudinal and vertical gusts. 


and the optimization had to progress by finding designs that passed this constraint. The 
pitch perturbation was also approaching the upper constraint limit in both the fixed and 
free static margin case, with the histogram indicating numerous designs violating the up- 
per constraint. An interesting behavior occurs in the cruise pitch perturbation response of 
Fig. 7.23(c). As the optimization progressed, the normalized constraint bounced back and 
forth between plus or minus 0.2. This behavior occurred due to some designs having the 
elevator response overshoot being more extreme than the initial deflection. As both ex- 
tremes were similar in magnitude, some designs had a positive value of the constraint while 
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Figure 7.23: SDynConFixSM and SDynConFreeSM configurations elevator response to pitch 
and airspeed perturbations. 


others had a negative value of the constraint, resulting in the back and forth behavior of 
Fig. 7.23(c). 

Examining the takeoff and maneuver condition static trim constraints, Fig. 7.24, the ele- 
vator effective angle of attack increased as the SDynConFreeSM optimization case progressed, 
whereas it remained nearly constant throughout the SDynConFixSM optimization. In both 
cases, the takeoff elevator constraint was active. The histograms of Figs. 7.24(b) and 7.24(d) 
show that nearly all the designs remained within the constraints with the SDynConFixSM 
case having more designs near the -11.4 degree constraint. As discussed previously, the 
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(a) SDynConFixSM normalized elevator deflections, (b) SDynConFixSM histogram of elevator deflections. 
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(c) SDynConFreeSM normalized elevator deflections. 



(d) SDynConFreeSM histogram of elevator deflec- 
tions. 


Figure 7.24: SDynConFixSM and SDynConFreeSM configurations trim deflection angles to 
the takeoff and maneuver flight conditions. 


maneuver condition resulted in positive lift on the horizontal stabilizer giving a positive 
effective angle of attack. However, neither the SDynConFixSM nor the SDynConFreeSM 
cases were constrained by the maneuver elevator constraint, with very few designs in the 
SDynConFreeSM case that resulted in angles greater than 12.2 degrees. 

Similar to the static trim constraint optimization cases described in Section 7.2, the one 
engine inoperative flight condition of Fig. 7.25 was a driving constraint on the vertical stabi- 
lizer sizing of the fixed and free static margin design optimizations. All aileron deflections, 
as before, were no greater than two degrees. Throughout the design optimization, shown in 


159 



Generation 



Deflection Angle, deg 


(a) SDynConFixSM normalized OEI deflections. (b) SDynConFixSM histogram of OEI deflections. 
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(c) SDynConFreeSM normalized OEI deflections. (d) SDynConFreeSM histogram of OEI deflections. 


Figure 7.25: SDynConFixSM and SDynConFreeSM configurations trim deflection angles in 
a one engine inoperative flight condition. 


the normalized deflection angles and histograms, the rudder deflection angles were within 
the constraint limits of plus or minus 20 degrees, but many designs were evaluated near the 
constraint. The best design rudder deflection angle increased throughout the optimization 
in the SDynConFreeSM optimization case reaching the constraint limit after 60 generations. 
Even though the total configuration geometry changed when the dynamic response con- 
straints were added, the one engine inoperative trim constraint was the driving constraint 
behind the sizing of the vertical stabilizers. 

Performing a multidisciplinary design optimization, including both static trim and dy- 
namic response constraints, produced feasible designs that included stability and control 
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considerations. Shown in Tables 7.12 and 7.13, none of the constraints were violated for ei- 
ther the SDynConFixSM or SDynConFreeSM optimal designs. The elevator effective angle 
of attack was nearing the negative constraint at takeoff while the rudder was having large 
deflection angles in the one engine inoperative condition. The lateral discrete gust placed 
a large demand on the aileron and rudder in the optimal SDynConFreeSM stall flight con- 
dition, but remained within the 20 degree deflection limit. Both the optimal designs easily 
met the roll residual constraint in both the cruise and stall flight conditions with the heavier 
roll weighting. The pitch residual performance was much improved from the static trim 
constraint optimization cases. 

Table 7.12: SDynConFixSM and SDynConFreeSM configurations static trim constraints. 


Condition 

Response Type 

SDynConFixSM 

SDynConFreeSM 

OEI 

Aileron 

-0.28 

0.28 


Rudder 

-16.62 

-19.62 

Takeoff 

Elevator 

-10.78 

-11.23 

Maneuver 

Elevator 

2.71 

-1.66 


Table 7.13: SDynConFixSM and SDynConFreeSM configurations dynamic response perfor- 
mance, 


Condition 

Response Type 

SDynConFixSM 

Cruise 

Stall 

SDynConFreeSM 

Cruise 

Stall 

Discrete Gust 

Aileron 

0.00 

8.36 

0.00 

-17.36 


Rudder 

-2.39 

-3.23 

-3.73 

17.59 


Elevator (Long) 

-2.04 

2.97 

-2.72 

1.19 


Elevator (Vert) 

-2.12 

3.82 

-2.87 

-1.31 

Turbulence 

RMS Pitch 

0.07 

0.72 

0.19 

0.96 


RMS Roll 

1.31 

0.99 

1.72 

1.24 


RMS Yaw 

0.39 

0.84 

0.45 

0.86 

Perturbation 

Aileron 

0.00 

12.06 

0.00 

12.45 


Rudder 

5.98 

5.32 

5.84 

5.46 


Elevator (Pitch) 

-3.00 

8.36 

3.31 

7.31 


Elevator (Speed) 

-2.55 

11.81 

-2.84 

12.11 

Residual 

Pitch 

-0.23 

-0.23 

-0.28 

-0.44 


Roll 

-0.53 

-0.47 

-0.56 

-0.48 


Speed 

0.00 

0.00 

0.00 

0.00 
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It was expected that the vertical gust would be sufficient to capture the largest deflec- 
tions placed on the elevator, and this appears to be the case. Table 7.13 shows that the 
longitudinal gust, for both the SDynConFixSM and SDynConFreeSM cases, placed a lower 
demand on the elevator in both the cruise and stall flight conditions. Although small in 
difference, this indicates the longitudinal gust could be excluded in the analysis, reducing 
both the computation time and number of design constraints. 

As the residual roll angle was consistently violated by the optimization cases not includ- 
ing the dynamic response constraints, and an additional roll angle weighting was added. Roll 
angle time response plots have been included in Fig. 7.26. The roll response was very similar 
for both configurations in the cruise and stall flight conditions. As indicated in Table 7.13, 
all cases had a roll residual angle within one degree of steady state five seconds after the 
initial perturbation. This can be seen in Fig. 7.26, as all the response lines were well above 
the intersection of the lower allowable constraint and the five second vertical dashed line. 
With the heavier weighting on the roll angle, the roll response has little to no overshoot in 
response to the perturbation, giving an excellent roll response. 
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(a) SDynConFixSM configuration cruise condition, (b) SDynConFreeSM configuration cruise condition. 



(c) SDynConFixSM configuration stall condition. (d) SDynConFreeSM configuration stall condition. 


Figure 7.26: SDynConFixSM and SDynConFreeSM configurations roll residual time re- 
sponses in the cruise and stall conditions. The horizontal dashed lines indicate the maximum 
allowed deviations after five seconds, highlighted with the vertical dashed line. 
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Chapter 8 

Summary and Conclusions 


A methodology for integrating stability and control into the conceptual design pro- 
cess, specifically focusing on multidisciplinary design optimization, was presented. Focus 
was placed on creating a controller that guaranteed closed-loop stability, providing system 
robustness during an optimization, while eliminating the necessity of the user being inti- 
mately involved with the controller design. This research broke away from the traditional 
military specifications on handling qualities, MIL-F-8785C and MIL-STD-1797A, where lin- 
ear dynamic modes must be identified, and requirements on mode natural frequencies and 
damping specified. Instead, this research used SAE-AS94900, which specified requirements 
on the system state response to perturbations, continuous turbulence, and discrete gusts. 
This eliminated the need to identify each linear mode during an optimization, which is es- 
pecially important for configurations, both with and without closed-loop control, where the 
traditional linear modes may not be present. Relying on the state response allowed for 
quantitative analysis of the system handling qualities, perfect for integration into a system 
optimization. The atmospheric disturbances and perturbations, along with static control sys- 
tem constraints, allowed for the geometric sizing of the system using actual system response 
requirements instead of the empirical tail volume coefficients often used during conceptual 
design. Key contributions of this research are summarized as: 

• Integrated stability and control into the conceptual design process in ModelCenter® 
using AVL, MATLAB®, and NASA’s Flight Optimization System (FLOPS) 

• Implemented a linear quadratic regulator (LQR) controller into the conceptual design 
process with gains selected using optimal control techniques 
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• Applied both static trim and dynamic response constraints to a multidisciplinary design 
optimization that negated the need to identify specific linear modes, focusing on state 
variable responses 

• Performed a multidisciplinary design optimization integrating mission performance us- 
ing FLOPS, aerodynamic analysis using AVL, and stability and control with a linear 
quadratic regulator stressed by atmospheric disturbances and perturbations 

Optimal control theory implementing a time-weighted linear quadratic regulator (LQR) 
was used for the closed-loop control design. The LQR methodology employed combined 
multiple techniques described in Ref. [4], but the complete derivation of the LQR technique 
used was described in Section 4.2, with additional information and derivations provided 
in Appendices B and C. The set of fully-coupled, linearized perturbation equations used 
were developed in conjunction with this research. Thrust terms and negligibly small terms 
were assumed to be zero, but greater knowledge of the system aerodynamic and thrust 
models will only improve the accuracy of the system dynamics. Increased accuracy in the 
information provided to the methodology will improve results as the technique is subject to 
the quality of information given to it. Sections 6.1 and 6.2 described the verification of the 
system dynamics and atmospheric disturbances using multiple publicly available resources 
including Refs. [15,29,30,78]. 

A Cessna 182T model was used as a test case throughout the methodology and analysis 
development. The geometry was taken from Ref. [15] and used to develop an AVL model. 
Techniques for modeling the fuselage and propeller effects were discussed, and the AVL 
model was validated using known stability data, presented in Refs. [15,34], A simple trade 
study was performed by varying the horizontal and vertical tail areas, which verified that the 
controller, the atmospheric disturbances, and the perturbations were working as expected. 

Described in Section 3.1, the D8.2b geometry was defined and used in the creation of 
OpenVSP and AVL models. Changes from the D8.1 geometry, reported in Ref. [10], were 
discussed. Experimental data for the D8.2b geometry was unavailable, so Cart3D, an inviscid 
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Euler solver, was used as validation of the D8 AVL model. From the aerodynamic validation 
study, modifications to the AVL model were made to better match the aerodynamic results 
of the Euler analysis. The modified AVL model was then used in a multidisciplinary design 
optimization. NASA’s Flight Optimization System (FLOPS) was used to size the D8 and 
estimate the weights, mission profile, and performance, and was integrated using ModelCen- 
ter with the flight dynamics and controls analysis and AVL. Phoenix Integration’s genetic 
algorithm, Darwin, was used in the system optimization, running five optimization cases, 
each with an increasing number of design variables and constraints. 

Overall, using fixed tail volume coefficients produced the expected design; the tail mo- 
ment arm was maximized to minimize the stabilizer area, resulting in reduced viscous drag 
and structural weight. This was accomplished by moving the wing apex to the forward most 
limit and sweeping the vertical stabilizer to the upper most limit, shifting the horizontal 
stabilizer to the most aft location. Adding the static trim constraints at takeoff, maneuver, 
and one engine inoperative flight conditions allowed the stabilizers to be sized using the 
static trim constraints instead of a fixed tail volume coefficient. This resulted in the wing 
apex shifting aft and the horizontal and vertical stabilizer areas decreasing dramatically, 
achieving a minimum total fuel burn and gross weight in the free static margin with static 
trim constraints optimization case, StaticConFreeSM. Takeoff trim and OEI constrained the 
elevator and rudder sizes, respectively. 

LIsing only the static constraints produced unreasonably small stabilizer areas that were 
increased with the addition of the dynamic constraints, placing a greater demand on the 
control system and stabilizer sizing. Nose-up trim at takeoff remained a constraining flight 
condition for the elevator, but the airspeed perturbation also became an active constraint. 
The elevator size was increased for the cases optimized with the dynamic constraints com- 
pared to the static trim constraints only. One engine inoperative remained the active sizing 
constraint on the rudder even with the dynamic constraints added. Takeoff and landing 
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field lengths were active constraints for all design optimizations except SDynConFreeSM, 
specifically sizing the wing area and static thrust. 

It was surprising that in the cases where static margin was allowed to be a free variable, 
the optimizer consistently chose to increase the static margin to take advantage of the natural 
stability benefits while decreasing the horizontal and vertical stabilizer areas. This indicated 
that the benefits from relaxed static stability, within the design space specified, were less 
than the weight and viscous drag benefits of downsizing the tails and shifting the neutral 
point/center of gravity, resulting in increased static margin. If center of gravity can be 
shifted as desired, the stabilizer areas should be decreased and static margin used to provide 
additional stability to take full advantage of the weight and viscous drag reductions. 

The residual pitch and roll angles after a system perturbation were active constraints in 
the design space, indicating that the design was seeking more authority from the controller. 
The weighting on the controller was selected to have an unbiased weighting on the state 
variable error; one degree of error was weighted the same as one foot per second of error. An 
additional weighting on the roll angle had to be added to improve the optimized solutions. 
It is hypothesized that additional performance gains could be achieved through adjustments 
to the LQR weighting matrices and time weighting. Adjusting the weighting matrices could 
allow a controller to have increased authority where needed, alleviating some of the active 
constraints. 

It was expected that the static trim constraints would have a greater impact on the 
horizontal and vertical stabilizer areas. In the one engine inoperative case, the engines 
were embedded in the fuselage, placing them near the center line. This resulted in a small 
moment arm in an engine-out case, requiring very little rudder authority. This allowed the 
twin- vertical stabilizers to be decreased far more than expected. It would be interesting to see 
how the OEI static trim constraint would affect the design with a podded or wing-mounted 
engine configuration. 
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Nose-up trim during the takeoff roll was far less constraining than anticipated on the 
optimization cases with only static trim constraints. The optimizer manipulated the geome- 
tries to place the neutral point and center of gravity in such a way as to have the wing lift 
counteract the nose-down pitching moment, reducing the required load to lift the nose wheel 
at the end of the takeoff roll. This was enabled by the center of gravity being dependent upon 
the neutral point placement and specified static margin. This flexibility in the placement 
is not usually feasible in conventional configurations but may be enabled in future concepts 
with large weight and density subsystems, such as batteries. 

Reviewing the disturbances used to stress the control system, it was seen that the 
RMS turbulence was never an active constraint in the design space. It was thought that 
the turbulence analysis would provide a unique constraint on the design space, eliminating 
designs that otherwise would have been acceptable. This did not occur. No cases evaluated 
ever had an RMS value over the constraint. It is possible the turbulence constraint could 
become active if a higher turbulence intensity was used, but this may become more of a 
structural problem than a stability and control problem. 

As anticipated, the crnise flight condition had little impact on the design space. When 
the dynamic constraints were not applied, the cruise condition roll residual was greater in the 
cruise condition than the stall condition. When the dynamic constraints were applied, both 
the stall and cruise residuals met the constraint and were close in value. It is believed that 
the stall residual constraint would capture the influence of the crnise roll residual constraint, 
eliminating the need to perform a dynamic analysis at the crnise condition. This would allow 
the addition of a new flight condition, perhaps a descent condition, that could add a new 
set of constraints on the design space without increasing the computation time. 

Reviewing the results, it is clear that going from fixed tail volume coefficients to static 
trim constraints, and then to static trim and dynamic constraints, drastically changed the 
optimal designs. Wing area and static thrust were nearly unchanged between the different 
optimization cases as they were driven by the takeoff and landing constraints, applied equally 
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to each optimization case. All other design variables changed dramatically as stability and 
control constraints were added to the design space. This shows the importance of including 
stability and control early in conceptual design, allowing it to influence the design space, 
thus taking advantage of the benefits it can provide instead of introducing penalties late in 
the design phase. The SDynConFreeSM case had reduced fuel burn compared to the baseline 
D8.2b, and the SDynConFreeSM configuration boasts an unblemished record of passing all 
constraints. It could be argued that the baseline design could pass the dynamic constraints 
with a controller redesign, as no static trim constraints were violated. However, the analysis 
performed here indicated that a horizontal stabilizer upload was required at the nose-up 
takeoff condition, even at the forward center of gravity position. This indicated that as the 
airplane accelerates down the runway and lift production increases, the natural tendency 
would be for an uncommanded pitch rotation counteracted only by a nose-down control 
input. As currently designed, this is perceived as a safety concern. 

Overall, including the dynamic constraints beyond the static trim constraints placed a 
great influence on the design space while showing that fixed tail volume coefficients and static 
trim constraints do not appropriately size the stabilizers. The dynamic constraints resulted 
in increased fuel burn compared to the StaticConFixSM and StaticConFreeSM cases, but the 
static trim only constraint cases produced infeasible designs when looking at stability and 
control considerations. The fuel burn for the SDynConFreeSM case was reduced compared 
to the baseline and fixed tail volume cases with additional reductions anticipated through 
adjustments in the LQR weighting matrices. Including stability and control in conceptual 
design only benefits the design, allowing for design space exploration in ways not possible 
when using the typical conceptual design disciplines. 
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Chapter 9 
Future Work 

Two studies that could be performed using the methodology as presented include an 
exploration of the system sensitivity to the system disturbances and an optimization study. 
The methodology presented used moderate disturbances in the discrete gust and continuous 
turbulence analysis, as specified in SAE-AS94900 [39]. The disturbance level increases as 
the probability of exceedance decreases, going all the way to an extreme disturbance. It 
would be interesting to see if the optimal designs shift with the introduction of higher 
disturbance intensities. As the disturbance probability decreases, the disturbance shifts 
to a more structural concern than a stability and control concern, assuming the disturbance 
does not cause an upset condition. If the configurations remain largely unchanged, the 
large disturbances should be considered mainly as a structural problem. However, if the 
configurations change drastically, considerations may have to be made to include these large 
disturbances in the configuration design or to resolve any upset potential in future detailed 
controller design. 

An optimization study would give greater confidence that the optimal solutions found 
were indeed the solutions that minimize fuel burn while meeting the applied constraints. 
Modifications to the population size, crossover, and mutation percentage would all have to 
be explored and any differences in results compared. This optimization study would require 
long computation times but would give greater confidence in the optimal solutions obtained. 
A current shortcoming is the system memory required for this type of optimization study, 
which results in halted optimizations from a combination of a large population size and a 
large number of generations. 
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Improvements to any analysis methodology can always be made by incorporating higher- 
order analysis methods, given the increased computation time is acceptable and the higher- 
order tools correlate to a higher hdelity solution. This research focused on maintaining a 
conceptual design level of configuration knowledge, catering each component to maintain 
that level of detail. Future capabilities that would increase the capability of the discussed 
methodology include: 

• Weight and balance calculations with higher-order mass moment of inertia calculations 

• Increase the number of flight conditions used to stress the configuration and flight 
control system placing additional constraints on the design space 

• Include elements of the performance output weighting matrix, W, and the control 
inputs weighting matrix, R , as design variables in the optimization 

The center of gravity was intentionally used as a dependent variable of neutral point 
location and static margin as part of the design space exploration in order to understand the 
system sensitivities. However, center of gravity is not a component that can be simply placed. 
The aircraft components and subsystems must be placed and configured to shift the center 
of gravity to the desirable position, making static margin dependent upon the configuration 
and subsystem layout. Static margin would become an optimization constraint instead of a 
design variable, restricting the allowable distance between the neutral point and center of 
gravity. This would require a weight and balance buildup needing either detailed placement 
of each component and subsystem or empirical data on each subsystem’s center of gravity. 

It has been emphasized that including stability and control in the conceptual design pro- 
cess is beneficial to the configuration design. A challenge of including flight dynamics early 
in the design is the detailed information required for accurate analysis, detail not typically 
available in conceptual design. Calculating highly accurate mass moments of inertia requires 
detailed component layouts and weight statements, much of which is unknown, undecided, 
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or fluidly changing. Implementing a method of calculating reasonably accurate mass mo- 
ments of inertia using conceptual design level details would be extremely beneficial to the 
methodology presented. Detailed design information would be unknown, so it is envisioned 
the calculations would use a semi-empirical methodology. Conventional and unconventional 
configurations, ranging from hybrid wing bodies to distributed propulsion concepts, would 
have to be included to give the greatest flexibility in exploring the design space. 

Discussed previously was the possibility of removing or replacing the dynamic analysis 
of the cruise flight condition, as it was not uniquely constraining on the design space. Any 
constraint violations in the cruise condition were also captured by the stall condition, making 
the cruise condition a potential unnecessary analysis. Adding a new flight condition to the 
analysis, such as a descent condition with flaps in the landing configuration, may impose new 
constraints on the design space. If computation time is not a concern, the new flight condition 
can simply be added to the current analysis with an additional set of constraints. This may 
provide a challenge to the optimization because currently a large number of constraints exist 
relative to a small number of design variables. In this case, there were 37 constraints to 9 
design variables (8 if static margin was fixed). The other option would be to replace the 
cruise condition with a new flight condition such as the descent flight condition. 

The cruise condition was only constraining on the design space in the pitch and roll 
residual constraints. It is hypothesized that the pitch and roll residual constraints could be 
alleviated through a modification of the LQR weighting matrices, placing greater weight on 
the pitch and roll states. The LQR weights were originally chosen as an unbiased weighting 
on each of the states as it was unknown what influence particular weights would have on 
each design. It is desired to add elements of the LQR weighting matrices as design vari- 
ables in the design optimization, giving greater controller design authority to the optimizer. 
It is envisioned that this expansion could provide additional benefits to the methodology 
presented, potentially improving the optimal designs and fuel burn reductions. 
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Even though minimization of analysis run time was a point of emphasis, further im- 
provements could be made to reduce the total system optimization time. The controller 
methodology and atmospheric disturbances were all coded using MATLAB®, which re- 
quires licensing and additional overhead in ModelCenter®. Using a compiled language, such 
as Fortran or C++, would allow for ModelCenter to simply call the executable without hav- 
ing to obtain a license, and performance improvements from a compiled versus interpreted 
language could be obtained. 

Replacing the vortex lattice aerodynamic analysis with a panel code would increase the 
order of the aerodynamic analysis, better capturing the actual geometry being modeled. A 
risk of using a panel code is the potential for increased computation time with the model 
being overly sensitive to panel size, especially as the configuration deviates from the baseline 
design during an optimization. Care would have to be taken to ensure the surface mesh 
generation was adequately automated to provide reasonable results (being better than the 
vortex lattice). This would be especially important for configurations that deviate from the 
traditional tube-and-wing concepts with conventional propulsion systems. 

Even though many configuration concepts are being explored for the future generation 
advanced air transport, concepts for the future 737 class air transport remain in the tube-and- 
wing configuration [10,83,84], These advanced concepts work well with the current method- 
ology as traditional control effectors were used, including ailerons, rudders, and an elevator. 
In these configurations the longitudinal and lateral/directional controls are mostly decou- 
pled, meaning the elevator is used primarily for longitudinal control and the ailerons and 
rudder for lateral/directional control. As currently implemented, the gain matrix contains 
all the cross-coupled terms in a 36-element matrix, resulting in a large number of variables 
minimized by the simplex algorithm. In the tube-and-wing concepts, the lateral/directional 
controls have negligible impact on the longitudinal controls. The gain matrix could possi- 
bly be reduced in size, removing the insignificant terms and, thus, reducing the number of 
variables in the simplex algorithm for minimizing the LQR performance index. This should 
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accelerate the LQR gain calculation convergence, saving computation time throughout the 
entire optimization. 

Decoupling the gain matrix would only be applicable for concepts with decoupled con- 
trols. For current military aircraft such as the F-22 and F-35, this would not work as their 
flight control systems actuate numerous effectors in response to a given input command, 
requiring cross-coupled control gains. Future air transport concepts, such as a blended wing 
body, will employ this control structure as well, requiring an expansion of the equations of 
motion used here. Non-unique control deflection solutions for redundant controls will also 
have to be addressed as researched by Garmendia in Ref. [85]. 

Distributed propulsion and boundary layer ingestion may require the equations of mo- 
tion to be modified to capture the unique aerodynamic-propulsion coupling provided by 
those technologies. Their effect on the concept stability control is not known, making his- 
torical data useless for conceptual design of stabilizing surfaces. The methodology presented 
here, with the appropriate system sensitivities captured in the equations of motion and sta- 
bility derivatives, would be perfect for assessing the stability and control characteristics of 
advanced concepts employing these technologies. The challenge would be in verifying and 
validating that the enhanced equations of motion accurately capture the complex coupling 
propulsion-airframe effects. 
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Appendix A 

Dynamic System Full Matrix Definitions 


The matrices of the fully coupled, perturbation equations are provided below. These 
matrices are used in Eqs. 4.1, 4.3, and 4.6. 
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The symbols used in Eqs. A.1-A.9 are defined as 
b = wingspan 

qoo = dynamic pressure 

a = steady-state angle of attack 

R = steady-state roll rate 

I' x , I' Y i I' z = stability axis mass 

moments of inertia 
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wing reference area 
aircraft mass 
gravity constant 
mean aerodynamic chord 
stability axis cross- 
product of inertia 
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Appendix B 

Derivation of the Standard PI State-feedback Constraint Equations 

The derivation of Eqs. 4.22-4.24 is given in this Appendix. Equations 4.20 and 4.21 
from Section 4.2.1 are repeated here for clarity. 

g = A t c P + PA c + Q + K t RK = 0 (B.l) 


Jtr = tr(PX) + tr(gS) 


(B.2) 


For the derivation of the state-feedback constraint equations some matrix calculus properties 
are required and are given as [4] 
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Taking the partial derivative of Eq. B.2 with respect to P gives 
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Expanding the last term on the right side of Eq. B.6 to get 
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which simplifies to 
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as given in Eq. 4.23. Taking the derivative of the Hamiltonian with respect to the matrix of 
Lagrange multipliers is straightforward resulting in the original constraint equation, g, given 
in Eq. B.9. 

fj 

Us^ds |tr (PA ' ) + tr {gS)] = 9 (B ’ 9) 

The final derivative of the Hamiltonian is taken with respect to the gain matrix, K 
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(B.10) 


where each term will be examined separately and are numbered accordingly. The term 
containing Q is independent of K and therefore has been neglected. Looking at the first 
term of Eq. B.10 and expanding for Hj 


tr (AjPS) = tr {A - BK) T PS = tr ( A T PS ) - tr (K t B t PS) (B.ll) 


The first term on the right hand side of Eq. B.ll is independent of the gain matrix goes to 
zero as the partial derivative is taken with respect to K . The second term on the right hand 
side of Eq. B.ll is a function of both K and K T requiring the use of the matrix property 
shown in Eq. B.5. 

A [_ (r (X T S T PS)] = [-tr (K t B t PS)] } T = -B t PS (B.12) 

Repeating the same methodology of Eqs. B.5 and B.12 on the second term of Eq. B.10 yields 
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Finally, the third term in Eq. B.10 is a function of both K and K T and so the product rule 
must be employed using the matrix calculus properties above to obtain 
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Summing the right hands side of Eqs. B.12, B.14, B.15, and B.16 and dividing by two gives 


0 = YM = RKS ~ BTps 


(B.17) 


which is the final constraint equation for the standard performance index linear quadratic 
regulator of Section 4.2.1, Eq. 4.24. 
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Appendix C 

Derivation of the Modified PI State-feedback Constraint Equations 


A cost function with time-dependent weighting is defined as 

J = 7: [ (t k e T e + u t Ru + z T Wz) dt 
2 Jo 


with the following: 
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Substituting Eq. C.2 into Eq. C.l gives 
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Splitting the cost function into two integrals 
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Using the method of integration by parts where u = t k , du = kt^ k ~^dt, dv = x T Px, and 
dv = x r Px and define — ^(x t Pqx) = x T Px, Eq. C.4 can be rewritten as [86] 
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For an asymptotically stable system this simplifies down to 


1 /*°° i roc 

J=- kt k ~ 1 x T P 0 xdt + - x t (K t RK + A T c H T WHA c )xdt (C.6) 

2 Jo 2 j 0 

The first iteration of the cost function in Eq. C.6 is subject to the constraint [4] 


^-(x T P Q x) = x T Px 

LLL 


(C.7) 


and using the definition x = A c x in Eq. C.7 gives the first of the nested Lyapunov equations 


LLjPn + P 0 A C + P — 0 


(C.8) 


This is the first constraint equation in a series of nested Lyapunov equations that were shown 
in Eq. 4.31. Continuing the process for another iteration, the method of integration by parts 
is again employed, using Eq. C.6, with u = kt^ k ~ l \ du — k (k — 1) t( k ~ 2 ^dt, dv = x T P 0 x , and 
v = — x T PiX, which defines — ^ (a; T P\x) = x T P 0 x. This gives the second iteration of the 
cost function as 


J = 



kt^ k ^ x T PiX 



(C.9) 


which, for an asymptotically stable system, simplifies to 


1 /»OG -| POO 

J=- k(k- 1) t^xTPrxdt +- x T ( K t RK + A t c H t WHA c ) xdt (C.10) 

2 Jo 2 Jq 

As in the previous iteration, Eq. C.10 is subject to the constraint 


- 1 


x T P 0 x 


(C.ll) 
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Substituting x = A c x into Eq. C.ll results in the second nested Lyapunov equation 


A]: Pi + P\A C + Pq — 0 


(C.12) 


The process of integrating the performance index by parts can be repeated until the following 
is obtained 


1 /»oo ^ ± 1 /*oo 

J = - Yl(k-n) x T Pk-ixdt +- x t (K t RK + A^H T WHA c )xdt (C.13) 

Jo n=0 Jo 


which can be simplified down to 


1 r°° 

J=- x T (k\P k - 1 + K t RK + A^H T WHA c )xdt (C.14) 

2 Jo 

Define 

- j^PkX = x T (k\Pk~i + I< t RK + A t c H t WHA c )x (C.15) 

When the left hand side of Eq. C.15 is substituted into Eq. C.14 and evaluated at the limits 
for an asymptotically stable system, the result is 


J = ^x T (0)P fc a;(0) 

Using x = A c x in Eq. C.15, the last Lyapunov equation is found to be 
A T c P k + P k Aj + k\P k _! + K t RK + A^H t WHA c = 0 


(C.16) 


(C.17) 


This is the last of the nested Lyapunov equations which are constraints to the solution to 
the modified performance index, Eq. C.16. The nested Lyapunov equations are summarized 
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below which matches Eq. 4.31 of Section 4.2.2. 


0 — go = AjPo + PqA c + P 
0 = 9i = Al^P\ + P\A C + Pq 

: (c.i8) 

0 = gk - i = A^P k _i + P k _\A c + P k -2 
0 = g k = AjP k + P k A c + fc!P fe _, + iL T AiL + A T C H T WHA C 

The Hamiltonian for the modified performance index is defined as [4] 

^ tr (P k . X) + htr (g 0 S 0 ) + ^ tr {giSi) H b^tr {g k -iS k _i) + ^ tr ( g k S k ) (C.19) 

where the derivative with respect to each independent variable P t , Si, and K is taken, where 
% — 0, 1, 2, . . . , k — 1, k. Taking the derivative of the Hamiltonian with respect to the matrix 
of Lagrange multipliers, So, Si, . . ., S k _ i, and S k results in 

dtf? 

- = 9i , i = 0,1,2, ... ,k (C.20) 

o Si 

which are the nested Lyapunov equations given in Eq. C.18. 

First, taking the partial derivative of the Hamiltonian with respect to Pq gives 

I f) 

0 = -QpP = 2 dP 0 = ^ tr <y9oS °^ + tr ( 9lSl ^ ( C - 21 ) 

where the terms independent of Pq have been removed for clarity. Expanding and simplifying 
Eq. C.21 results in 

H c £ 0 + S 0 Aj; + S 1 = 0 (C.22) 
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The process repeats following the same pattern as Eq. C.22 until the k — 1 term is reached 


giving 


(9 W 1 

0 = ttw — = A c Sk-i + Sk- i^4j + k\S k 


dP k - 1 

Finally, taking the derivative of the Hamiltonian with respect to P k gives 


(C.23) 


834? 

m 


o =^r = A c S k + S k A^ + X 


(C.24) 


The nested Lyapunov equations of Lagrange multipliers is summarized as 


0 — A c S k + S ’fcHj + X 
0 = A c Sk~ i + + k\Sk 


(C.25) 


0 = A c Si + .SiHj + S 2 
0 = H C S 0 + S 0 AJ + 5! 


Taking the derivative of the Hamiltonian with respect to the gain matrix, K , is best 
shown by breaking down Eq. C.19 into the individual terms. The second term on the right 
hand side of Eq. C.19 is the first term dependent on K . The first partial derivative gives 

\~^K tr (k9qS ^ = \~0 k (" 4 ' rjPo5 ’°) + tr (Wo) + tr (PS'o)] (C.26) 


The last term on the right hand side of Eq. C.26 is independent of K and is therefore zero. 
The closed-loop state matrix is dependent on K and must be expanded. The expanded right 
hand side, neglecting the last term in Eq. C.26 produces 


ld_ 
2d K 


tr 


(A — BK) t P 0 S 0 


+ tr 


P 0 ( A - BK) So 


(C.27) 
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Using the matrix properties given in Eqs. B.3-B.5, the solution to Eq. C.27 is 

~t r (g 0 S 0 ) = -B T P 0 S 0 (C.28) 

This process can be repeated for all terms of the form \tr ( gjSi ) in the Hamiltonian, for 
i — 0, 1, , k — 1, with the solution being 

~tr ( 9i Si) = -B T P l S i , i = 0,l,...,k-l (C.29) 

Looking at the last term of the Hamiltonian, the k th term, and neglecting terms independent 
of K gives 


lm tr{9kSk) = 

1 d 


2M tT ^ ( A c PkSk "> + tr ( PkA cSk) + tr ( K Tr KS k) + tr ( AjH T WHA c S k )] (C.30) 


The first two terms of Eq. C.30 will give —B^P^Sk, which is of the same form as Eq. C.29. 
Taking the derivative of the third term is quite simple and results in 


1 d 


[tr ( K T RKS k )] = \ ( RKS k + RKS k ) = RKS k 


(C.31) 


The last term of the Hamiltonian needs to be expanded and factored in order to take the 
partial derivative with respect to K . First, substituting for the closed-loop state matrix gives 

[Jk [ tr (AlH T WHA c S k )} = 1 -P [tr [(yl T - K T B T ) H r WH ( A - BI< ) Sj } (C.32) 
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Expanding the right hand side of Eq. C.32 results in four terms to look at individually as 


~4r-tr f A T H T WHAS k -A T H T WHBKS k - 
2 oh \ 


© 


K B H W H AS k + KBH WHBKS k 


(C.33) 


The first term in Eq. C.33 is independent of K and therefore goes to zero. Taking the 
derivative of the second and third term results in 


l ~tr (■ -A T H T WHBKS k ) = ~B T H T WHAS k (C.34) 

l ~tr ( -K T B T H T WHAS k ) = ~B T H T WHAS k (C.35) 


Using the chain rule in the fourth term, the derivative is found be 


-J^tr ( K T B T H T WHBKS k ) = B T H T WHBKS k 

2 oh v ' 


(C.36) 


Summing all the components gives the full solution 


~dK 


RKS k - B t ( P 0 S 0 + P 1 S 1 + --- + P k -uS k -i + P k S k ) 

+ B T H T WHBKS k - B T H T WHAS k (C.37) 


When using a simplex algorithm for the minimization of the performance index, only 
Eqs. C.16 and C.18 need to be used. However, when a gradient-based minimization routine 
is used, Eqs. C.16, C.18, C.25, and C.37 need to be used, where Eq. C.37 is used as the 
gradient. 
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